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The  work  covered  by  this  repoit  was  done  cm  the  Sy'stcm  Dynamics  Branch, 
Aeronautical  Research  Laboratory,  under  Proa-rt  7060,  "Flight  Dynamics 
Research  and  Analysis  Facility".  Mr.  Paul.  .v_  Xoskcr  is  Project  Engineer. 

This  study  is  part  of  a  continuing  program  to  n^ermine  optimum  methods  of 
simulation  and  analysis  of  the  dynamics  of  ai~  weapon  systems.  The  general 
subject  of  quaternions  as  applied  to  coordinate  conversions  has  been  under 
investigation  for  approximately  two  years,  thcmch  the  bulk  of  the  work  reported 
here  was  accomplished  during  the  last, six  moir.cn s  of  1957. 

The  author  wishes  to  express  his  appreciation  to  Mr.  Robert  T.  Harnett  and 
others  of  the  Analog  Computation  Branch  of  tint-  -Aeronautical  Research  Laboratory 
for  assistance  in  the  analog  simulation  portiair.  nl  the  study. 
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ABSTRACT 


The  theory  of  the  four -parameter  method  is  developed  with  specific 
application  to  coordinate  conversion  in  aircraft  simulations.  'I  his  method 
is  compared  with  the  direction  cosine  method  both  in  a  theoretical  error 
analysis  and  in  an  example  simulation  on  an  analog  computer.  It  is  shown 
that  the  quaternion  method  is  no  more  sensitive  to  multiplier  errors  than 
is  the  direction  cosine  method,  and  it  requires  nearly  iO  per  cent  less 
computing  equipment.  In  addition,  the  multiplier  bandpass  requirement 
in  the  four -parameter  method  is  only  half  as  severe  as  for  direction  cosines. 
By  every  important  criterion,  the  quaternion  method  is  no  worse  than,  and 
in  most  cases,  better  than  the  direction  cosine  method. 
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SECTION  I 


INTRODUCTION  j 

The  problem  of  motion  of  a  rigid  body  and  the  associated  one  of  coordinate 
conversion  are  very  old  ones  in  the  field  of  classical  dynamics.  Significant 
results,  dating  from  the  time  of  Euler  (1776)  through  the  introduction  and 
application  of  matrix  methods  by  Cayley  and  Klein  and  others  in  the  last  half 
of  the  nineteenth  century,  brought  the  matter  to  such  a  satisfactory  state  that 
no  significantly  new  methods  or  approaches  have  been  found  necessary.  The 
development  of  modern  computing  machinery  makes  necessary  a  re-examination 
of  the  various  methods  from  the  standpoint  of  their  utility  in  computational 
devices.  It  is  not  necessarily  true  that  methods  which  have  proven  their  con¬ 
venience  in  the  largely  analytical  manipulations  of  classical  mechanics  should 
prove  to  be  best  adapted  for  numerical  or  a  lalog  computation.  Quaternions  fell  ( 

into  disuse  among  physicists  about  the  turn  of  the  present  centurv  because  matrix  i 

and  vector  methods  had  proved  more  useful  in  the  types  of  investigations  then 
being  conducted.  The  purpose  of  the  present  paper  is  to  show  that  the  quaternion  , 

appro  a  .:h  to  coordinate  transformation  does  offer  veal  advantages  in  the  analog 
simulation  of  rigid  body  motion.  In  recent  times  Deschamps  and  Sudduth*  have 
suggested  an  application  lot  digital  compulation,  and  Backus**  has  proposed  them 
for  analog  simulation,  but  in  general  quaternions  are  little  known  among  those 
engaged  in  simulation  of  aircraft  motions,  ! 

The  coordinate  conversion  problem  in  aircraft  and  missile  simulation  is 
different  at  least  in  emphasis  from  that  of  classical  dynamics.,  It  might  be  well  j 

to  state  the  problem  which  is  of  interest  and  to  which  the  methods  explained  later 
will  be  applied,  A  missile  or  aircraft  tnay  be  considered  as  a  moving  coordinate  '■ 

system.  Various  vectors  must  be  transformed  into  this  coordinate  system  or  out 

♦Deschamps,  G.  A.  and  W.  B,  Sudduth,  Federal  Telecommunications  Laboratories, 
Nutley,.  New  Jersey,  Case  <16-10707,  November  1955. 

•v* backus,  George,  Rigid  Body  Equations  -  Euler  Paramete?  s.  Technical  Note  6, 
Advisory  Board  on  Simulation,  University  of  Chicago,  November  1951. 
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of  it  into  semi:  inertial  system.  Integrating  the  equations  o.r  motion  of  the  air¬ 
frame  can  be  made  to  yield  the  three  components  of  the  coordinate  system's 
angular  velocity  vector.  From  the  X,  Y  and  Z  components  (P,  Q,  R)  of  this 
vector  in  the  moving  system,  it  is  desired  to  keep  track  of  the  orientation  of 
the  coordinate  system  in  such  a  way  that  vectors  may  be  transformed  in  either 
direction.  This  means  an  integration  of  angular  rate  to  determine  angular 
position. 

Fundamental  to  this  procedure  is  a  consideration  of  how  the  orientation  of 
the  coordinate  system  is  to  be  specified.  During  the  history  of  the  subject, 
various  methods  of  doing  this  have  been  put  forward.  Alt  the  most  useful  ones 
fail  into  three  categories:  Euler  angles,  quaternions,  and  direction  cosines. 

Of  these,  the  first  and  last  are  probably  the  most  familiar  to  modern  readers. 

In  the  Euler  angle  method,  the  orientation  is  expressed  as  the  result  of  three 
rotations  about  each  of  three  axes,  the  rotations  being  made  in  a  specific 
sequence.  The  physical  interpretation  of  a  quateinion  is  a  rotation  through  some 
angle  about  some  specific  fixed  axis.  The  nine  direction  cosines  are  simply  the 
cosines  of  the  angles  between  each  of  the  axes  in  the  moving  system  with  each 
of  the  axes  of  the  fixed  system.  Principal  attention  here  will  be  given  to  the 
quaternion,  or  four -parameter  system.  It  was  first  introduced  by  Euler  in  1  776, 
a;  a  rt-suU  of  spherical  trigonometry  considei  ations.  The  elegant  quaternion 
formulation  was  invented  by  Hamilton  in  1  843  as  a  new  kind  of  algebraic  formal¬ 
ism.  A  matrix  formulation  was  devised  by  Klein  for  use  in  gyroscopic  problems 
and,  in  this  formulation,  is  usually  known  as  the  Caylcy-Klein  parameters.  Each 
cf  these  three  different  approaches  to  the  four-parameter  system  has  its  own 
advantages.  It  has  been  decided  to  present  at  least  an  outline  of  ail  three  here. 
There  are  two  reasons  for  this:  first,  there  are  some  propositions  which  arc  more 
easily  shown  by  one  development;  second,  it  seems  probable  that  when  the  reader 
is  offered  a  choice  oi  method,  he  will  reach  an  under  standing  sooner  it  he  con 
select  the  method  most  nearly  consonant  with  his  own  backgiound. 

It  will  become  apparent  that  this  subject  presents  something  of  an  expositional 
problem.  In  order  to  reach  the  desired  ends,  it  has  been  deemed  to  assume  that 
the  reader  has  a  knowledge  of  matrix  methods,  especially  as  applied  to  coordinate 
conversion  in  three-dimensional  space.  As  a  compromise,  a  brief  introduction 
to  the  subject  is  given  in  Appendix  A,  though  a  more  satisfactory  treatment  is 
given  by  Goldstein*.  In  this  i  eport  the  term  "quaternion"  lias  been  used  to 


d.-lein,  Herbert,  Classical  Mechanics,  Addison- vVesley  Press,  (.-air.bri  dSe-‘» 
Mass.,  1?50. 
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represent  the  four -parameter  method  in  general.  In  other  cases,  it  is  necessary 
to  use  the  word  to  distinguish  Hamilton's  development  from  the  others.  It  is  hoped 
that  confusion  may  be  kept  to  a  minimum. 

There  arc  many  different  techniques  used  in  present-day  aircraft  simulations  to 
solve  the  coordinate  conversion  problem.  The  technique  is  usually  adapted  to  the 
special  requirements  of  the  problem  at  hand.  If  most  of  the  rotation  takes  place 
about  one  axis,  or  if  only  the  gravity  vector  is  to  be  handled,  or  if  the  airframe's 
rotation  is  otherwise  restricted,  valuable  simplifications  may  be  effected  in  the 
analog  equipment  required  to  represent  the  conversion.  It  is  not  the  present 
purpose,  however,  to  investigate  all  these  possibilities.  Consideration  will  be 
given  only  to  the  most  general  and  unrestricted  case:  that  of  several  complete 
revolutions  about  any  cr  all  axes.  This  immediately  excludes  the  Euler  angles 
because  of  the  singular  point.  The  advantages  of  Euler  angles  are  such,  and  their 
popularity  is  so  pervasive,  however,  as  to  warrant  keeping  them  in  mind.  Accordingly, 
Appendix  B  gives  a  brief  outline  of  the  Euler  Angle  system  most  commonly  used  in 
aircraft  work,  and  at  appropriate  points,  comparisons  will  be  made  of  them  with 
quaternions  and  direction  cosines.  In  making  such  comparisons,  that  form  of  Euler 
angle  instrumentation  whose  capabilities  most  nearly  equal  those  of  the  quaternion 
scheme  will  be  assumed.  This  form  has  been  discussed  at  some  length  by  Howe* 
and  his  figures  and  results  will  be  used  for  comparison.  In  Howe's  method,  the 
extent  and  direction  of  rotation  is  unrestricted  except  for  the  inevitable  singular 
orientation,  and  he  shows  that  even  this  leads  to  less  practical  difficulties  than  one 
might  expect. 

It  is  valuable  to  keep  the  Euler  angles  in  mind,  but  the  quaternion  method  must 
really  stand  or  fall  or.  its  comparison  with  direction  cosines.  It  has  in  common 
with  direction  cosines  the  capability  of  handling  completely  unrestricted  rotations. 
Accordingly,  considerable  attention  has  been  devoted  to  the  direction  cosine  method 
in  this  report.  Both  a  theoretical  error  analysis  and  a  simulation  program  were  done 
for  the  cosines  in  order  tc  provide  the  most  complete  possible  basis  of  comparison. 
They  have  been  dune  before,  bet  it  is  difficult  to  compare  results  obtained  by 
dilferent  investigators  or.  difft  ’-ent  computing  equipment.  An  attempt  was  made  here 
to  keep  the  conditions  as  nearly  comparable  as  possible.  Of  all  the  material  con¬ 
tained  herein,  no  originality  is  claimed  except  for  the  quaternion  error  analysis 
and  simulation.  Even  here,  no  ne.v  techniques  were  used,  with  the  possible 
exception  of  tbe  method  of  handling  multiplier  errors.  It  was  felt  necessary, 
however,  to  include  the  remaining  material  in  order  to  introduce  and  place  in 
context  this  probably  unfamiliar  subject. 

♦Howe,  R.  M.  and  E.  G.  Gilbcii,  A  New  Resolving  Method  for  Analog  Computers, 
WADC  Technical  Note  55-468,  January  1  956. 
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SECTION  II 


•bt 


THE  EULER  PARAMETERS 


The  earliest  formulation  of  the  four -parameter  system  was  given  by  Euler 

1  776,  though  the  oldest  treatment  generally  available  today  is  probably  that 
of  Whittaker*.  It  is  an  essentially  geometrical  development,  but  will  not  be 
presented  as  such  here.  The  principal  results  may  be  demonstrated  with  much 
less  labor  by  use  of  matrices. 

Central  to  the  development  of  these  parameters,  and  indeed  to  the  four- 
parametei  methods  in  general,  is  the  proposition  known  as  Euler's  theorem, 
which  may  be  stated  as  follows:  any  real  rotation  may  be  expressed  as  a 
rotation  through  some  angle,  about  some  fixed  axis.  In  other  words,  regard¬ 
less  of  what  the  rotation  history  of  a  body  is,  once  it  reaches  some  orientation, 
that  orientation  may  be  specified  in  terms  of  a  rotation  through  some  angle  (which 
can  be  determined)  about  some  fixed  axis. 

The  truth  of  this  proposition  is  not  intuitively  obvious,  but  in  any  case,  it  must 
be  shown.  Consider  a  transformation  matrix  (A).  No  restrictions  are  pur.  on  (A) 
other  than  those  which  exist  for  all  orthogonal  transformation  matrices  (see 
Appendix  A),  Another  way  of  stating  Euler's  theorem  is  to  say  that  for  every 
matrix  (A)  there  exists  some  vector  R  whose  components  are  the  same  before 
and  after  application  of  (A);  in  other  words,  there  must  be  rorr.e  R  such  that 

(A)R  =  R  (1) 

for  any  (A),  If  the  components  of  R  are  designated  X,  Y  and  Z,  the  elements 
of  (A)  by  ,  then  Equation  (1)  may  be  written 


all 

a12 

aI3\ 

x\ 

\ 

X' 

i 

;  a21 

a22  ' 

a23 

r  - 

y 

\> 

a32 

a  / 
a33  / 

/ 

V 

\  z 

If  this  matrix  equation  is  expanded  in  components,  a  set  of  linear  homogenous 


♦Whittaker,  E.  T.  Analytical  Dynamics,  Fourth  Edition,  Dover  Publications, 

N.  Y.,  I'M 4. 
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equations  results: 


(a,  j  -  1)X  +  a12Y  +  a,3Z  =  0  , 

a21X  +l*h  '  ‘>Y  +  a23Z  =  °*  (3) 

a^jX  t  a.,Y  +  (a^  -  1)Z  =  0. 

A  necessary  and  sufficient  condition  for  existence  of  a  non-trivial  solution  is 
that  the  determinant  of  coefficients  be  zero.  Therefore,  it  is  necessary  to  show 

that 


11  -1 

*12 

a13 

a2 1 

*22  '  1 

a23 

a31 

a32 

a33  -  1 

This  may  easily  be  done  making  use  of  the  properties  of  an  orthogonal  trans¬ 
formation  matrix  developed  in  Appendix  A.  If  the  above  equation  is  expanded, 

(al  la22a33+a12a23a31+a21a32a13“a31a13a22'a21a12a33'a32a23ai  1 

+(ai  1  -a22?-3  3+a2  3a3z)+(a2^'al  la3  3+al  3a31^a3  3"a  1 1*22**21*12^  °  ^ 


The  first  term  vanishes  in  consequence  of  the  fact  that  the  determinant  of  the 
transformation  matrix  must  equal  unity  (Eiquation  (156)  ),  and  the  last  three  terms 
vanish  from  the  orthogonality  conditions  of  Equation  (162).  Thus,  it  is  proved 
that  Equation  (4)  is  an  identity  tor  any  orthogonal  (A)  and  that  there  exists  some 
vector  R  which  is  unchanged  by  the  transformation.  This  proves  Eulerls  theorem. 

Since  it  has  been  shown  that  it  is  possible  to  express  any  rotation  as  a  single 
rotation  about  some  axis,  it  is  possible  to  make  use  of  the  equivalent  rotation  to 
specify  orientation.  Consider  two  coordinate  systems  X  Y  Z  and  X*  Y’21.  The 
XYZ  system  is  assumed  to  be  fixed  in  inertial  space,  and  X’Y'Z'  is  moving  in 
some  arbitrary  manner,  though  both  coordinate  systems  have  the  same  origin. 
Assume  that  initially  the  two  systems  are  coincident.  Then  the  X'Y'Z.'  system 
is  rotated  through  an  angle  p  about  an  axis  which  makes  angles  n,  p,  Y  with  the 
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X,  \,  /.  axt-s  res pectivcly.  Ji  will  bo  noted  that  tins  axis  of  rotation  makes  the 

s.-.iiu;  .ingles  a,  p,  \  with  the  X1,  Y1,  Z*  axes  also.  It  is  now  necessary  to  express 

the  t  i  .-.nsformution  matrix  in  terms  of  the  quantities  pi,  q,  p  and 

In  order  to  do  this,  use  is  made  of  an  additional  coordinate  system,  X  Y  Z  , 

1  r  r  r 

ithn  h  is  fixed  in  the  XY  Z  system.  The  axis  lies  along  the  axis  of  rotation, 

.ind  the  Y'  axis  is  restricted  to  the  XY  plane,  'ibis  would  give  rise  to  difficulty 

it  the  Z  axis  is  the  axis  of  rotation,  but  in  that  case,  the  Y  axis  could  be  confined 

to  the  XZ  plane  or  the  YZ  plane,  and  the  final  result  would  be  unaltered.  At  any 

rx.te,  with  the  choice  indicated,  tlie  Y  axis  is  always  perpendicular  to  the  Z  axis.. 

Now  the  rotation  through  the  angle  p  is  a  rotation  through  p  about  the  X^  axis,  so 

the  rotation  is  a  very  simple  one  in  the  system.  Accordingly,  the  rotr'ian 

o!  ’he  X'Y'Z1  system  th  ’ough  the  angle  p.  may  be  viewed  as  the  result  of  three 

rotations:  (1)  rotation  of  the  X'Y'Z'  system  into  coincidence  with  the  X^Y^Z^ 

system;  (2)  rotation  through  the  angle  p  abc  t  the  axis;  (3)  the  reverse  of 

(1)  to  restore  the  original  separation  of  the  X'Y'Z'  and  X  Y^Z^  systems.  The 

matrix  for  each  of  these  transformations  will  be  developed,  and  then  the  three 

may  be  multiplied  together  to  express  the  total  ti  ansfor  mation. 

Nirsi,  the  transformation  into  the  X  Y  Z  system  will  be  considered,  a,  p, 

r  r  r 

and  v  are  the  angles  between  the  now  X  axis  and  the  fixed  X,  Y  and  Z  axes. 
Thus,  it  is  seen  from  Eq  nation  (1Z5)  that  ajj,  aJ2  and  are  immediately  fixed. 
One  other  cosine  may  be  established.  Recall  that  the  Y^  axis  is  perpendicular  to 
me  Z  mxis.  Ibis  means  that  a^  =  0*  Thus  the  matrix  of  the  first  rotation  is 
partially  <■  st ablishcd, 


/cos  a 

cos  P 

II 

a2  2 

\a31 

a3Z 

Applying  the  orthogonality 

conditions. 

«.  laments  are 

cos  y 

0 

*33 

it  is  possible  to  deduce  that  the  other 


cos  a 


cos  p 


cos  v 


(A) 


r  cos  p  escY 
+  cos  a  cot  v 


i  cos  a  esc  \  0 

+  cos  p  cot  y  ±  sin  y 


(V 


The  ambiguities  m  sign  may  be  resolved  by  making  use  of 
’’  matrix  above  must  reduce  to  the  identity  matrix  when 


the  requirement  that 
a  becomes  zero. 
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The  result  is 


cos  a  cos  p  cos  y 

-cos  p  esc  v  cos  a  esc  y  0 

-cos  a  cot  y  -cos  p  cot  y  sin  y 


The  second  rotation,  through  the  angle  p,  about  the  axi3  is  simply 


(R)  = 


0 

cos  p. 
-sin  p. 


0 

sin  p. 
cos  p. 


The  last  of  the  three  rotations  is  the  inverse  of  (A)  or  (A)”  .  Thus,  the  general 
transformation  is  the  result  of  all  three,  called  (B).  It  is  given  by 

(B)  »  (A)-’(R)(A).  (10) 

This  is  a  similarity  transformation,  and,  among  other  things,  the  spur  (sum  of 
the  diagonal  elements)  of  a  matrix  is  invariant  under  a. similarity  transformation, 

i.  e. ,  '■  " .  .  .  .  \  ' 

>  bj  j  +  *  ^33  , =  1+2  cos  p,  (11) 

so  the  angle’ of  rotation  may  be  obtained  directly  from  the  diagonal  elements  of 
the  transformation  matrix.  Carrying  out  the  operations  of  Equation  (10)  gives 


1  -2  sin^ 

u  .  2 
sin  a 

2(sin^  ^  cos  a  cos  p 

2(cos  a  cos  y  sin^  ^ 

. 

-sin cos  cos  y) 

-sin  cos  ^  cos  p)  ’ 

2(sin^4jf  cos  a  cos  p 

1  -2  sin^  sin^p 

2 

2(sin  ^-cos  p  cos  y 

-sin-^-  cos  cos  y) 

-sin  cos-^-  cos  a) 

^  2(cos  a  cos  y  sin^-^-  2(sin^-£  cos  p  cos  y  1-2  sin^-^-sin^  y 


. +  sin-jp-cos  ^  cofr(V)  -sin4£  cos -jy  cos  a) 

If  the  following  substitutions  are  made,  ^ 

|  =  cos  a  sin-^,  tj  =  cos  P  sin-t£>  §  =  cos  y  sin-tjr,  x  =  cos-^,  ^  (13) 


WADC  TR  58-17 


the  matrix  of  i  1 )  bit  omo.- 


B  - 


r  V  -i*  •  X 

Hi  n  - 4  x) 

Hi  4  *■  n  x) 


<£(  6*1  + ;  x) 

.  I  .2,2 

*4  +  r\  -4  +  X 

2<l4  -ix) 


2(4  4  -  ix) 

2(*li+  4  x) 

,2  2  .2  ^  2 

-4  -*)  +4  +x 


(U) 


These  four  quantities  are  called  the  Euler  parameters.  It  may  be  seen  from 
their  definitions  that  they  obey  the  relationship 


.2,  21.21  2  , 

4  +  n  +4  +  x  »ii 


(15) 


so  they  are  not  all  independent.  Also,  none  may  lie  outside  the  range  *1. 

If  the  quantities  p,  a,  (3,  and  y  are  known,  it  is  a  simple  matter  to  compute 
the  Euler  parameters  and/or  the  transformation  matrix  by  the  method  given  above. 
If,  on  the  other  hand,  the  transformation  matrix  is  given,  it  is  also  possible  to 
solve  for  the  four  parc^neters,  though  difficulties  arise.  A  consideration  of  these 
difficulties  will  shed  further  light  on  the  nature  of  the  Euler  parameters.  To  begin 
with,  it  should  be  stated  that  the  quantities  p,  a,  p,  and  y  cannot  be  uniquely 
determined  from  the  transformation  matrix.  The  reason  for  this  is  that  even  though 
rotation  through  a  certain  angle,  about  a  certain  axis  will  produce  a  definite  unam¬ 
biguous  orientation,  the  reverse  ie  not  true.  If  the  orientation  is  given,  there  are 
four  separate  ways  in  which  it  could  have  been  obtained  by  rotation  about  a  fixed 
axis.  Possibly  a-'  example  will  help  to  clarify  this.  Assume  that  the  rotation  being 
considered  is  a  rotation  through  an  angle  of  +  30°  about  the  +  X  axis.  There  are 

three  other  ways  to  get  to  the  same  position:  (1)  a  rotation  through  -  30°  about  the 

- - - - -  />  .  . 

-  X  axis;  (2)  a  rotation  through  -  330°  about  the  +  X  axis;  (3)  a  rotation  through 
^  about  the  -  X  axis.  A  further  illustratlon-of  the  possibilities  is  given  in 

the  table  following. 


Case  1 
Case  L 
Case  3 
Case  4 


X 

4 

n 

4 

+  cos-jr 

+  cos  a  sin^- 

+  cos  P  sin 

+  cos  y  sin-jj- 

+  cos 

c 

(-  cos  a)(-  sin  ■£) 
+  cos  a  {-  sin 4“) 

(-  cos  a)  sin  £ 

(-  cos  P)(-  sin  •£) 

+  cos  p  (-  sin  •£) 
cos  P)  sin  ^ 

(-  cos  y)(  -  sinJjr) 
+  cos  y  (-  sin  -£) 
(-  cos  y)  ,sin  £ 

-  cos£ 
u 

-  cos£ 

The  first  two  cases  lead  to  the  same  Euler  parameters,  and  the  last  two  lead  to 
a  different  set  which  are  the  negative  of  the  first.  All  four  sets  lead  to  the  same 
transformation  matrix. 
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The  relationship  between  Euler  parameters  and  direction  cosines  may 
be  derived  by  equaling  terms  ir.  Equation  rli).  The  result  is 

-  If  an  +  a2Z  +  a33, 

4§  -  1  +  aj  j  -  a22  -  a33,  (16) 

2 

4t|  =  1  -  an  +  a2Z  -  a33. 

4C,  =  1  -  aj  3  -  a,2  4  a3J. 

These  eqtiations  determine  the  Euler  parameters  except  for  sign.  The  sign  must 
be  gotten  in  another  way.  From  comparison  of  terms  in  the  matrix  it  is  possible 
to  show  that 

a31  ~  a?  3  =  4*t1* 

a12  -  *21  =  4X4.  (17) 

a23  -  a32  =  4*^ 

Thus,  if  x  ts  assumed  to  be  always  positive,  the  signs  of  the  others  may  be 
deduced  from  Equations  (17)  unless  x  =  0*  This  is  the  special  case  of  a  180° 
rotation.  There  is  an  additional  (tiubiguity  hei  e  oecause  tne  eirection  of  ihe 
axis  of  rotation  and  the  direction  of  the  rot  ition  arc  completely  unrelated. 

Either  a  positive  or  a  negative  rotation  about  either  the  positive  or  negative 
axis  will  give  the  same  result.  For  this  special  case,  another  means  would 
have  to  be  devised  for  defining  the  signs,  but  it  hardly  seems  worthwhile  to  go 
into  it  here.  It  is  net  expected  that  this  will  lead  to  any  practical  difficulties. 
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SECTION  III 


THE  CAYLEY -KLEIN  PARAMETERS 

In  this  development  of  the  four-parameter  system,  it  is  found  that  a 
2x2  complex  matrix  may  be  used  to  represent  a  real  rotation,  rather  than  a 
3x3  real  matrix.  Consider  such  a  matrix  (H), 


(H)  = 


The  requirement  is  placed  on  this  matrix  that  it  be  unitary,  that  is  to  say  the 
product  of  (111  and  its  adjoint  must  yield  the  unit  matrix.  The  adjoint  is  the 
complex  conjugate  of  the  transposed  matrix.  In  addition,  it  is  required  that 
the  determinant  of  the  matrix  (II)  have  the  value  +1.  The  unitary  condition 
allows  ±1  for  the  determinant,  so  this  is  an  additional  requirement.  The 
unitary  condition  may  be  written  as 


1 1 

h12 

21 

h22 

hll*  hzi* 
h12*  h22* 


hll  b12 
hZl  h22 


1  0 
(  0  1 


Expending  and  equating  components  gives 


11 

*h 
nI  l 

+ 

b21 

*hi 

=  1, 

11 

*k12 

+ 

h21 

‘hi 

-  0, 

12 

*hn 

+ 

h22 

*h21 

=  o. 

12 

*hl  2 

+ 

b2Z 

‘hi 

=  1. 

The  second  and  third  equations  arc  the  same,  being  merely  complex  conjugates 
of  each  other.  The  fir3t  and  fourth  equations  h-ave  no  imaginary  component, 
whereas  the  second  (or  third)  has  both  real  and  imaginary  parts.  Therefore,  the 
three  independent  equations  contain  four  conditions.  These,  together  with  the  deter 
minant  requirement  that  hjjh^  -  =  m ak e  it  possible  to  determine  Certain 

relationships  among  the  four  quantities  ^nin*  ^  maY  be  shown  that  =  hjj* 


and  li 


,j  -  -  hj,*  so  (lie  matrix  may  bo  written  as 
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(H) 


hll  n12 


-h  *  h  * 
n 1 2  n\l 


(21) 


The  quantities  h^,  hj^,  h ^  are  usually  referred  to  as  the  Cayley-Klein 
parameters.  It  will  be  noted  that,  they  are  complex  numbers.  While  it  is 
convenient  to  use  them  as  such  ’  ’  analytical  operations  (and  this  is  the 
purpose  for  which  Klein  developed  then:)  a  physical  computer  must  treat 
complex  numbers  in  terms  of  their  real  and  imaginary  parts.  Therefore,, 
it  is  convenient  to  introduce  fou  *  other  quantities  defined  as  follows: 


hn 

=  ei  +  ie 

hl  2 

=  e3  +  ie 

(22) 


where  the  e's  are  all  real  numbers,  and  i  is  the  square  roof  of  -1.  Using 
these  definitions,  the  matr  ix  (H)  may  be  written  as 


el  +  ie2 


e3  +  ie4 


(H)  = 

h  "  ie2  i 

Now  consider  another  complex  matrix  (P),  which  has  the  form 


(2  3) 


(P) 


x  -  ly 


x  +  iy 


(24) 


where  x,  y  and  z  arc  real  numbers.  It  will  be  noted  that  the  matrix  (P)  is 
equal  to  its  own  adjoint,  and  thus  is  said  to  be  self-adjoint  or  Hermitian.  Nov/ 
consider  a  transformation  of  (P)  of  the  form 

(P)'  =  <H)(P)(H)+  (25) 

where  (H)+  designates  the  adjoint  of  (hi).  Since  (H)  is  unitary,  (H)*  -•  (H)”*, 
so  equation  (2  5)  is 

(P)'  =  (H)(P)(H)"1  .  (26) 


This  is  a  similarity  transformation.  It  is  shown  in  Appendix  A  that  the  deter¬ 
minant  of  a  matrix  is  invariant  under  a  similarity  transformation.  It  can  also 
be  shown  that  tire  Hermitian  property  and  the  spur  are  both  invariant  under  a 
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til^u 


similarity  transformation.  Therefore,  the  transformed  matijx  (P)'  must 
have  the  form 


x'  -  iy ' 


<pv-  |  • 

\x'  +  iy'  -z'  ‘  / 


The  fact  that  the  determinant  cf  (P)  must  equal  the  determinant  of  (P)'  gives 

x2  +  y2  +  z2  =  x'2  +  y’2  +  z’2.  (28) 

It  x,  y  and  a  are  viewed  as  components  of  a  vector,  then  Equation  (28)  is  the 
requirement  that  the  length  of  the  vector  remain  unch.-mged.  Equation  (26)  may 
be:  written 


z1  x'-iy1 
x'-Hy'  ~z' 


:1+iez  e3hie4\  j  z  x-iv  \  I  e}  -ie2  -e3-ie4 
-e3t-ie4  ej  -i  ..  )  ^x+iy-a  j  \  e3-ie4 


If  the  operations  of  Equation  (29)  a  *e  carried  out,  it  is  found  that 

x!  =  -  e3"  +  c4“)x  -  2'e^  +  03e4)y  x  2(e2«4  -  e^)55  < 

>■'  =  2(c3e4  -  CjU^Jx  +  (ej2  -  cz''  e32  -  e4”)y  +  2(8^  +  e^jjz  ,  (30) 

2  2  2  2' 

z'  =  2 ( e  j  e 3  »•  e^e4)>;  +  2(e?e'3  ••  e,e4)y  +  (m  +  -  c3  -  ')2  . 

Those  equations  represent  a  linear  transformation  between  the  components  of  x,  y 
and  z,  and  the  components  oi  x1  y1  and  z1.  The  matrix  for  this  transformation 


222,2  -y  t  .  \  \  \ 

el  -c2  ~e3  -a4  ?(elc2+<;3e4)  2(,i2e4'e!e3)  \ 


(A)  =  2(e3c4-e?e?<) 


2(clcd^ 


2  222  ,  > 
;1  "e2  '"1-3  -e4  2'e2e3+e4el' 


7.(eie3-e,e4) 


2  2  2  2 
el  ie2  'c3  ” e4 


It  may  be  shown  directly  that  tills  matrix  satisfies  tne  orthogonality  conditions, 
but  it  is  proved  alto  from  Equation  (23).  Equation  (31)  shows  that  the  nine 
direction  cosines  may  be  expressed  in  terms  of  llit:  four  e's.  If  Equations  (22) 
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vittil  into  L^utiUunb  (*.0j  it  is  found  that 


2  l  Z  Z  . 

'l  +  ei  +  *3  +  *4  *  *•  <-,-J 

and  therefore,  only  three  of  the  er*  haditipendenl.  The  identity  o!  the*  c 

four  quantities  with  the  Euler  parameter*  is  obvious.  Comparison  of  Equations 
(31 )  and  ( 1 4)  gives 

el  ~  Xj  3  l.  ”  n.  %  =  4-  (33) 

An  equivalence  ha*  been  indicated  between  the  roai  (3x3)  matrix  (A)  and  the 
complex  (2x2)  matrix  (If)-  it  may  b#  shown  that  this  correspondence  goes 
further.  Consider  the  real  transforms  tit  n 

.  £  «  (B)?. _ .... 

and  let  the  associated  unitary  complex  matrix  be  IH)^,  ao  that 

(P)'  =  (H)1(R)(H)J'  . 

Nov/  consider  a  second  transformation  (A)  with  associated  (H)2» 

r"  =  <A)r'  . 

(p)"  r  (H)2(P)‘(H>2  , 

Substituting  (34)  and  (35)  Into  (36)  gives 

r"  =  (C)T. 

(”)"  =  -  <»7) 

Therefore,  if  (A)(B)  =  (C)  and  (H)2(H>,  =  (H)y  the  above  equations  become 

*’•  r  (C)r. 

(P)”  =  {H)5(P)(H)3  ,  (38) 

_ _ _  ♦  Vi  a  t  mu'tipHcafion  of  two  real  1x3  ratrlce#  correspond*  to  multiplication 

of  the  two  associated  2x2  complex  matrices  in  tie  twit  order.  Two  types  of  quan¬ 
tities  which  correspond  in  this  manner  ere  said  to  be  iBvmorphic. 


(34) 

(35) 

(36) 
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It  is  *1.0  po.fible  to  view  ttiis  proceee  o.»  two  succsosive  rotation,  in  terms 
toe  e1.  themselves.  Con. id. r  on*  rotation  defined  by  *v  e^,  Cj  and  e4- 
After  this,  another  rotation  is  performed  which  io  deucribed  by  t,j',  e^',  *3' 
»ade4',  There  ie  *om«  set  of  e'.  called  ej",  e^",  #j..#  u^.  wMch  de.cribc. 
the  final  orientation  after  the  two  rotation*.  This  combined  »»t  may  bo  found  by 
multiplying  ih*  (H)  matrices  of  the  two  rotation*  In  the  correct.  sequence,  The 
equation  in 


e3"HV\  /•l’«V 

’03"*iG*"  '  I  «rs 


Expanding  this  equation  anti  equating  component.  give. 


®l"  '  *l'®l  ~  *1*1  “  93  *3  -  «4'®4  * 

V  =  eZel’  ♦  *1*1  *  °3'“4  -  *4  *3  , 

•3"  r  °i'*s  ‘  °z’H  +  -3'«i  ♦  %'•*  •  (40> 

. .  V  “  *z‘«3  +  V«4  +  V*1  -  *l'*Z  ’ 

By  uae  of  these  equation.,  successive  transformations  may  ho  handled  In  terms  of 
the  e's  direct'y 

This  technique  may  be  used  to  determine  the  relationship  between  the  e'o  and 
the  Euler  angles  given  in  Appendix  B.  The  (H)  matrix  corresponding  to  each  of 
the  Euler  angle  rotations  may  be  determined,  and  tho  three  may  be  multiplied  in 
tho  correct  order  to  synthesize  the  complete  transformation.  Consider  first  the 
<H)  matrix  corresponding  1<j  tits  first  Euler  a&gls,  given  as  <!•  in  Appendix  Bc 
F-om  Equation  (179)  it  ia  seen  that  the  transformation  equations  are 


x’  V  x  cq«  4  •>  y  *iis  4-, 

y'  >;  -*  sin  +  +  y  co*  <j«r  (41) 

s’  =■  *, 

/liquating  coefficient,  of  these  equation.*  with  like  coefficient,  in  Equations  (50) 
give,  the  nine  relations 
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i 

4 


L.  .  U.. 


i  2  2  2,2  . 

cos  iji  =  -e,  -e3  +e4  ,  -  sin  4-  =  2(°3e .j-Cj  e^). 


0  :  ^(elc3^c2C4^' 


sin  ij/  =  2(eie24e3e4), 
0  —  2(e  ,e^  -ej  c^). 


2  2‘  3  eA  U  ^  2(c2e3-eie4),  (42) 

,  2a  2  2  2 
1  =  C1  +e2  'e3  -e4  \ 


cos  41  ~  -e^  +c*3  ~^4 

0  ~  ^(e2C3+e4ej)r 


Those  equations  cannot  all  be  satisfied  unless  e,  =  e.  -  0.  If  this  is  true,  then 
2  2  .  .  ,  J  4  2  2 


cos  4>.=  el  -  e2  . 


sin  4<  -  Zc^c^, 


el  '*2  =>> 


Solving  these  equations  for  e^  and  gives 


*|>  .  4* 

T  *  e2  =  slnT  * 

so  the  (H)  matrix  corresponding  to  the  41  rotation  is 


e ^  =  cos  -v-  , 


(43) 

(44) 


I  it)  ,  4* 

1  cos  -y  +  i  sin  -y 


\ 


4> 


cos  ~  i  sin 


I  4 

I  e  2 


e  2  0 


i  ^ 
0  e"1? 


(45) 


By  an  exactly  similar  process,  it  may  by  shown  that  the  other  two  matrices  arc 


(46) 


Therefore,  che  entire  tr ansforination,  which  is  the  result  of  all  three  rotations, 
is 


/  0 
fcos-7 
- 

•  0  \ 
siny  | 

i 

fcos$ 

isinf  ’ 

1  ■  9 

e 

.  (Hg  = 

i .  .  ± 

± 

! 

\  ~sm2 

cosz/ 

l  ism 

cos  f 

(H)  = 


+  ie? 


«  -A  -f  ’ 

\  ~3+--4 


e3  +  ie4 


e  -  -  ic- 

l  c  s 


(»)+(H)e(H)+  . 


(47) 


Carrying  out  the  indicated  multiplications,  and  equating  components  gives 


£ 

2 


* 

el 

i. 

=  cos  y 

0 

COS  J 

<b 

cos  f 

+ 

sin  | 

.  0 
sln  T 

e2 

.  ‘P 
=  Sln  7 

G 

COS  1 

cos  2 

- 

cosi 

0 

Sin  ry 

Cm 

c3 

•-  cos 

*  ^ 

.  e 
Sln2 

cos 

+ 

•  i 
sln  2 

G 

c0s  z 

e4 

4< 

=  COSj 

G 

cosT 

sin 

- 

stn  y 

.  G 
sm7 

(48) 


2 

± 

2 
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SECTION  IV 


QUATERNIONS 

The  most  brilliant  formulation  of  the  four-parameter  method  was  made  by 
Hamilton  in  1S43.  He  developed  a  new  type  of  entity  called  a  "quaternion".  It 
is  composed  of  four  parts, 

q  =  S  f  ia  +  jb  +  kc,  (49) 

where  S,  a,  b  and  c  are  real  numbers,  and  the  indices  i,  j  and  k  are  defined 
by  the  following  rules; 

i2  =  -  1 ,  ij  =  -  ji  =  k, 

j2  =  -  1,  jk  =  -  kj  =  i,  (50) 

k2  =  -  1,  ki  =  -  ik  =  j. 

The  conjugate  of  the  q  ’aternion  q  is 

q"'  =  S  -  ia  -  jb  -  kc.  (51) 

Using  the  laws  for  the  indices  quoted  above,  it  may  be  easily  shown  that 

qq*  =  q*q  =  S2  +  a"  +  b2  r  c2,  (52) 

which  is  called  the  length  or  norm  of  the  quaternion.  If  this  norm  is  unity,  then 
a  special  form  of  quaternion  results,  a  versor.  It  is  possible  to  make  use  of 
these  to  describe  a  coordinate  transformation.  The  quantity  S  is  called  the  real 
or  scalar  part  of  the  quaternion,  and  ia  J-  jb  +  kc  is  called  the  imaginary  or 
vector  part.  Now  assume  we  have  a  quaternion  whose  scalar  part  is  zero.  We 
call  this  a  vector  of  components  X,  Y  and  Z, 

V  =  iX  +  jY  +  kZ.  (53) 

Let  us  examine  the  operation 

q*Vq  =  V'  (54) 

where  q  is  a  versor.  So  far  there  is  no  particular  reason  to  expect  that  V 
will  be  a  vector,  but  this  turns  out  to  be  the  case.  Equation  (54)  may  be  written 

(S  -  ia  -  jb  -  kc)  (iX  +  jY  +  kZ)  (S  +  ia  -5-  jb  +  kc)  =  V'.  (55)  -- 

Wher  this  equation  is  expanded  making  use  of  the  rules  for  indices,  the  result  is 

V'  '  -  c2]  +  Y  [2cs+2ab]  +  Z[2ac-2sb]} 

L  x !_  <iao  -  2cs]  +  Y  f  f  S2  -  a2  -  b2  -  c2]  +  Z  [  2as  +  2cb]  }  (56) 

{  x  [  2  b  -  ^  r]  +  Y  [2bc  -  2sa]  +  Z  [s2  -  a2  -  b2  +  c2]  }. 

-C  TP 
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This,  is  simply  a  coordinate  transformation  whose  I  ransfor  jnalion  matrice  is 


2  2,2  2. 

s  4  a  -  b  -  c 


2(ab  -  cs) 
2(ac  4  sb) 


2(c.  s  4  ab) 

2  2  ,  .  2  2 
s  -  a  4-  d  -  c 

2 (be  -  sa) 


2(ac  -  sb) 
2(as  4  cb) 


2  .  2  2 
a  -  b  4  c 


The  correlations  with  matrices  derived  in  the  two  preceding  sections  ar; 
evidently 


ei  =  X  =  s. 


eL-l-~  c* 


e3  =  rj  =  b. 


«4  =  t  ~  a- 


The  matter  of  two  successive  rotations  may  be  handled  quite  easily.  Assume 
that  first  we  transform  a  vector  with  the  versor  q^. 


q^Vqj  =V*. 


Next  we  apply  the  versor  q^. 


V”  =  q2*V'qz  =  q^+Vq^.  (60) 

We  now  define  a  new  vector  q^q^  =  q  —  and  wish  to  find  the  relationship  between 
q3  and  q2*qi  4',  \V e  define  It  may  be  seen  that 

3  qzq4,  (61) 

ar.d  since  is  a  versor,  q^q^*  =  !•  Therefore,  Equation  (6l)  reduces  to 


^i  =  W 


Now  we  apply  on  the  left. 


.  .  <Uq =  <llV4  =  1  _=  q3q4  * 

so  that  must  equal  the  conjugate  of  q^.  This  means  that 

'V  =  q3*Vq3  . 

Now  observe  that  the  equation  q3  -  q^q^  may  be  written 


S3  -  ia3  —  jb3  "  he 3  5  ,Sj  4  b  jbj  ^  hCj)  (S^  4-  ia^  4  jb^  4  hc^). 
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Expanding  this  equation  and  equating  components  gives 

S3  --  S,S2  -  a2  “  blb2  '  clc2  * 

a3  =  S 1  a2  +  S2a]  +  blC2  "  Clb2  ’ 

(66) 

b3  “  Slb2  '  alC2  +  blS2  +  Cla2  ’ 

c3  =  S1C2  +  alb2  “  bla2  +  C1S2  * 

These  equations  are  identical  with  Equations  (40)  which  were  developed  in  the 
same  connection  by  use  of  the  Cayley-Klein  parameters.  Thus,  the  quaternion 
method  leads  to  the  same  result  as  the  preceeding  developments. 


I 
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SECTION  V 


1NFINIT  ESKIA  L  TRANSFORMATIONS  AND  RATE  OF  ROTATION 


The  preceding  sections  have  dealt  with  the  four-parameter  method  of 
specifying  the  orientation  of  a  coordinate  system.  As  was  stated  in  Section  I, 
however,  the  primary  interest  is  in  determining  the  orientation  from  the  rate 
of  rotation  through  a  process  of  integration.  Accordingly,  it  is  necessary  to 
relate  the  rates  of  change  ol  the  four  parameters  to  the  rates  oi  rotation  of 
the  axis  system. 

It  was  shewn  in  Section  III  that  an  orthogonal  transformation  may  be 
represented  by  a  complex  matrix  having  certain  properties.  It  is  now  of 
interest  to  investigate  this  matrix  when  an  infinitcsmal  rotation  is  performed. 
Let  us  assume  that  this  infinitesmal  rotation  consists  of  a  rotation  through  the 
angle  Ap  about  a  line  which  tnakes  angles  ol  a,  p  and  y  with  the  X,  Y  and 
Z  axes  respectively.  Recall  that  the  matrix  (H)  may  be  expressed 


(H)  = 


el  +  ie2  c3  +  ie4  \ 


c ,  Tie,  e.  -  1  e  , 

3  4  12 


Applying  the  geometrical  interpretation  of  the  e's  gives 


<H)  = 


cos-^y  Ti  cos  y  sin-^- 

k>  t- 


-  cos  p  sin  ty  T  i  cos  a  sin  iy 


cos  p  sin-t  +  i  cos  a  sin  -t 


cos  ^  ~  i  cos  y  si n  ^ 


From  this,  it  is  possible  to  see  that  the  infinitesmal  rotation  may  be  represented 

by 


IT  i  ^  cos  Y 


P  T  i  ^  cos  a 


cos  P  4  i  cos  a 
.  Ap 

l-i  cos  y 


since  cos  ~  1,  sin  ^  ~  . 


It  is  expected  that  any  matrix  representing  an  infinitesmal  rotation  will  differ 
only  slightly  from  the  identity  matrix.  This  is  true  of  the  above  matrix,  and 
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■his  may  be  shown  more  clearly  by  writing  it  as  follows: 


1  cos  y 


-  cos  p  4-  i  cos  a 


cos  {3  +  i  cos  a 


-l  cos  y 


=  (I)  -(O. 


Mow  assume  that  this  infinitesmal  rotation  takes  place  during  a  small  time  interval 
At.  If  (H)  is  the  matrix  at  the  beginning  of  the  interval,  and  if  (H)'  is  the  matrix 
at  the  end  of  the  interval,  then  the  time  derivative  of  (ll)  may  be  written  as 


_d  _  lirr 

dt  (H'  ”  At-*-C 


lim  (H)'  -  (H) 
=  At -Hi  At 


The  final  matrix  (H)'  may  also  be  viewed  as  the  result  of  two  rotations,  first  (H) 
and  then  (H)^  .  In  other  /ords,  (H)'  =  (H)f  (H).  Putting  this  into  the  above 


equation  gives 


i-  (">  -  A™0  #  (") 


Since  (H)  is  not  affected  by  the  time  inclement,  the  limiting  process  refers  only 
to  the  quantity 


(e  )  _  1  Au  i  cos  y 
At  2  At  -  cos  p  +  i  cos  a 


cos  p  t  i  cos  a 
-  i  cos  y 


I  :i  the  limit,  the  quantity  -^j;is  simply  the  scalar  magnitude  of  the  angular 
velocity  vector.  If  P,  Q  -  and  R  are  the  components  of  this  velocity  vector 
ab  ng  the  X,  Y  and  Z  axes,  then  evidently  cos  c  =  P,  cos  y  -  R, 
cos  p  =  Q,  so  that 


/  iR 

lim  (e )  _  1  I 
At  —  0  At  =  T  \  -  Q  +  iP 


Q  +  iP 


Therefore,  from  Equation  (72), 


ahH>4 


iR 

Q  *-  \? 


Q  +  iP 

-  iR 


is  a:  m  poss- 


,  a  straightforward  limiting  process,- that  the  time 
also  a  matrix  whose  elements  are  the  time  derivatives 
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of  the  elements  of  the  original  matrix.  Therefore, 


•  .  • 

el  +  le2 

e3  +  ie4  1 

1 

~  2 

f  iR 

Q  +  iP  | 

el  +ic2 

e3  +  ie4 

-;3  +  U4 

*i  -ihzj 

|  -Q  +  iP 

-mj 

~e3  +  ic4 

C1  '  ie2 

Expanding  and  equating  like  components  gives 


(76) 


Zel  =  "  *4P  ‘  e3°  "  e2R* 

2e,  =  -  e,P  +  e.Q  +  e.R, 

2  3  4  i  (77) 

2e3  =  +  e^P  +  ejQ  -  e4R, 

2e4  =  +  ejP  -  e2Q  +  e3R. 

These  are  the  equations  which  would  be  used  to  compute  the  four  parameters  in 
an  actual  simulation.  Now  if  we  multiply  Equation  (76)  on  the  right  by  the  adjoint 
of  (H)  the  result  is 


el  +ie2 


e3  +  ie4 


el  "  ie2 


e3  -  ie4 


“e3  +  ie4 


el  *  ie2 


e3  ‘  ie4 


el  +  ie2 


.Again  expanding  and  equating  components  gives 


r 

\  -Q  -  iP 


P  =  2  {-  e4ej  -  e3e2  +  e2e3  +  e^), 

Q  =  2  (-  e3ej  +  e4e2  +  e^  -  e2*e4)f  (79) 

R  =  2  (-  e2^  +  ej^  -  e4;3  +  e^). 

Thus,  if  the  four  parameters  and  their  rates  of  change  are  known,  the  angular 
velocity  may  be  computed. 
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SECTION  VI 


THEORETICAL  ERROR  ANALYSES 


In  the  preceeding  sections,  the  fundamental  theory  of  the  quaternion 
method  has  been  presented.  Before  proceeding  to  an  application  of  the 
method,  it  is  of  interest  to  study  theoretically  the  errors  to  be  expected. 

Not  only  will  this  give  a  prediction  of  the  results  to  be  obtained  in  the 
simulation,  but  it  will  give  a  better  understanding  of  just  how  the  equations 
must  be  instrumented  in  order  to  achieve  the  maximum  accuracy  of  which 
the  computer  equipment  is  capable. 

As  was  mentioned  earlier,  both  quaternion  and  direction  cosines  will  be 
simulated,  so  errors  for  both  were  analyzed  on  much  the  same  basis.  It  is 
felt  that  this  is  an  important  part  of  the  demonstration,  because  without  a 
theoretical  error  compari  on,  any  differences  found  in  the  simulation  would 
be  subject  to  the  question  o  computer  malfunction.  If  simulator  results  and 
theoretical  error  analyses  agree  with  each  other,  the  degree  of  confidence  in 
the  comparison  will  be  much  higher.  Theoretical  error  analysis  is  but  little 
used  by  analog  computer  operators,  especially  in  non-linear  problems  such  as 
this.  It  turns  out,  however,  that  both  quaternions  and  direction  cosines  lend 
themselves  readily  to  an  analysis  of  errors  and  the  results  obtained  agree  with 
observations. 


A .  Direction  Cosine  Method 


The  fundamental  equations  to  be  used  in  generating  the  direction  cosine 
t  ransfo  rmation  are  given  in  Appendix  A,  There  are,  however,  many  possible 
variations  which  will  be  discussed  briefly.  Possibly  the  most  straightforward 
way  would  he  to  solve  the  nine  simultaneous  equations  and  thus  generate  all  nine 
of  the  direction  cosines  by  integration.  As  the  solution  progresses,  however, 


it  is  inevitable  that  errors  will  accumulate.  Some  of  these 


v-iJ  VV  V  _  .  I 


.  o  i  r>  wiu  t/c  wi 


such  nature  as  to  cause  the  orthogonality  conditions  (Equations  130)  not  to  be 
satisfied  after  a  time  if,  indeed,  they  were  satisfied  initially.  This  may  be 
thought  of  as  a  departure  of  the  three  axes  of  the  moving  system  from  mutual 
orthogonality  and  distortion  of  the  unit  le-yth  of  these  axes.  Some  of  the  errors 
rising  ;n  the  solution  will  r-t  contribute  to  this,  and  these  may  be  thought  of  as 
mguiai  irif<-  o,  the  coordinate  system  will  drift  as  a  whole,  and 


V  3!  -i 
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in  addition,  the  unit  vectors  will  change  their  relative  orientations  and  lengths. 
These  last  two  difficulties  may  be  eliminated,  as  will  be  shown  presently,  but 
it  should  be  understood  that  this  is  of  but  little  value  unless  some  way  can  be 
found  for  making  the  drift  of  the  system  as  a  whole  tolerable.  In  an  aircraft 
simulation  where  a  coordinate  conversion  is  used,  there  generally  exists  some 
feedback  which  will  eliminate  long-term  drift  in  the  coordinate  system.  As  will 
be  shown  later,  the  drift  can  be  reduced  to  where  it  is  much  less  than  those  drift- 
producing  elements  in  the  physical  system  being  simulated,  such  as  airframe 
misalignment,  gyro  drift  and  amplifier  noise.  If  the  errors  due  to  rotational  drift 
cannot  be  corrected,  there  is  not  much  additional  penalty  in  accepting  the  errors 
due  to  non -orthogonality.  In  any  case,  it  is  advisable  to  determine  in  advance  how 
much  drift  can  be  allowed  in  the  given  application,  and  to  design  the  coordinate 
conversion  to  meet  the  requirements,  using  the  techniques  developed  in  this  section. 

The  possibility  of  correcting  orthogonality  errors  was  first  suggested  by 
Corbett*.  Possibly  a  description  of  the  corrections  in  physical  terms  will  be  the 
most  instructive.  It  may  be  seen  from  the  material  presented  in  Appendix  A, 
that  a  physical  interpretation  may  be  placed  on  the  rows  and  columns  of  the  trans¬ 
formation  matrix  (A).  The  elements  of  the  first  row,  for  instance,  may  be  consid¬ 
ered  as  the  three  components  of  the  unit  vector  i'  along  the  three  unprimed  (fixed) 
axes.  Similarly  the  elements  of  the  first  column  may  be  viewed  as  the  three  com¬ 
ponents  of  the  unit  vector  i  along  the  three  primed  (moving)  axes.  Both  i^  ^,-lc 
and  i',  j'  ,  k'  coordinate  vectors  are  orthonormal  sets.  These  facts  may  be 
written  as 

t*r=i  T* f  =  o  T*  •T*  =  i  T*  =  o 

M  =  1  f  •?  =  0  y  =  i  j1  *  p  =  o  (80) 

k  "k  =  1  k  *T  =  0  lT*  ’"k*  =  1  k'  *T'  =  o 

These  are  vector  equations,  and  may,  therefore,  be  expanded  in  any  coordinate 
system.  Expanding  the  first  six  in  the  primed  system  and  the  last  six  in  the  un¬ 
primed  system-gives  _  _ 


♦Unfortunately,  this  work  has  not  been  generally  available.  The  first  published 
document  is  WADC  Technical  Report  57-425  Stabilization  of  Computer  Circuits. 
November  1957. 
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**• 


2  .  2 
all  +  a21 

X  2 

+  a31  = 

l. 

2  .  2 
al 2  +  a22 

,  2 
+  a32  = 

1. 

2‘  2 
a13  +  a23 

.  2 
+  a33  = 

1, 

alla12  +  a2 1  a22 

+  a31a32 

=  o. 

a 1 1 a 1 3  +  a21a23 

+  a31a33 

=  o. 

al 2al 3  +  a22a23 

+  a32a33 

=  0, 

2  .  2 
all  +a12 

+  a132  = 

1. 

2  2 
a21  +  a22 

+  a232  = 

1. 

2  .  2 
a31  +  a32 

+  a333  = 

1» 

alla21  +  a12a22 

h  3a23 

=  0, 

alla31  +  a12a32 

+  aj  3a33 

=  0, 

a21a31  +  a22a32 

+  a23a33 

=  0. 

(81) 


The  first  six  of  these  will  be  recognized  as  the  orthogonality  conditions  given  in 
Appendix  A.  The  last  six  may  be  seen  to  be  the  same  conditions  for  the  matrix 
(A.)-* .  All  twelve  equations  must  be  satisfied  by  any  real  orthogonal  transform¬ 
ation  matrix. 

The  general  procedure  is  to  compute  the  direction  cosines  by  integration,  take 
the  computed- quantities  and  perform  the  operations  of  some  of  the  Equations  (81). 

If  the  equations  (81)  are  not  satisfied,  the  error  is  used  to  modify  the  direction 
cosines  until  they  are.  It  appears  to  be  necessary  to  compute  at  least  six  of  the 
direction  cosines  by  integration.  Several  schemes  were  tried  for  getting  by  with 
less  than  six,  but  none  of  those  e  stable.  Assume  for  instance  that  the  three 
omp'  .ients  ^ •'  ...cd  system,  and  the  three  components  of  j  in  the 

r>*-*  '  imputed.  The  three  components  of  k  may  then  be  computed 


;  7R 
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without  integrations  from  the  relationship  that  k  =  ixj.  If  six  direction  cosines 
are  computed  by  integration,  then  three  of  Equations  (81)  must  be  used  to 
eliminate  redundancy.  If  nme  <~csinc-s  are  integrated,  then  six  auxiliary  equation 
must  be  used.  Both  possibilities  were  proposed  by  Corbett*,  and  an  improved 
version  of  the  second  by  Howe**.  The  first  alternative  (called  the  two-vector 
method)  was  the  one  chosen  for  use  here.  It  was  selected  because  it  uses  less 
simulation  equipment,  though  the  advantage  over  Howe's  three -vector  version 
is  not  great.  Both  require  36  multiplications,  though  the  two -vector  method 
requires  only  six  integrators  rather  than  nine. 

The  two-vector  method  may  be  described  as  follows.  The  three  components 
of  i  in  the  primed  system  are  computed  by  integration  through  use  of  the  k^, 
and  a^j  differential  equations.  Then  a  normalizing  circuit  is  added  to  keep 
the  length  of  the  vector  unity.  In  addition,  the  components  of  j  in  the  primed 
system  are  computed.  These  are  a^,  azZ  anc*  a32*  ^  correction  is  added  to 

keep  this  vector's  length  equal  to  unity;  another  correction  is  added  to  keep  this 

vector  normal  to  the  first.  Finally,  the  components  of  k  are  computed  from 

— ►  —  ►  ■  ► 

the  equation  k  =  i  x  j.  The  complete  set  of  equations  to  be  solved  is 

all  =  a21R  '  a31°  ‘  klall  ^  ‘  1 

a21  =  a31P  "  allR  “  kla21  ^  "  1  1)» 

a21  =  allQ  '  a21P  ‘  kla31  ^  ~  i'i^ 

a12  =  a22R  “  a32^  ”  kla12  ^  "  J  j)  ~  k2al  *  j)» 

(82) 

a22  =  a32P  "  a12R  ‘  kla22  ^  j  ^  '  k2a21^ 

a32  =  a12Q  '  a22P  "  kla32  ^  j  ^  ‘  k2a31^ 

al 3  =  a2ia32  *  a22a3T 


a23  =  a12a31  ‘  al la32* 


a33  =  alla22  ‘  a12a21* 


♦Corbett,  op.  cit. 

**Howe,  R.  M.  Coordinate  Systems  and  Methods  of  Coordinate  Transformations 
for  Dimensional  Flight  Equations  Proceedings  of  the  First  Flight  Simul; 
tion  Symposium,  November  56,  WSPG  Special  Report  9,  White  Sands 
Proving  Grounc. 
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where 


t.t  Z  Z  .  Z 
j  J  =  a,2  +  aZ2  +  a32  , 


2  2  .  2 

1  1  =  all  +  '21  +  a31 

i  *j  -  ana12  +  a21a22  +  a31a32. 

It  will  be  observed  that  there  are  six  differential  equations,  and  three  algebraic 
equations,  so  that  six  integrators  will  be  required.  The  kj  terms  are  for 
correction  of  the  length  of  the  vectors  and  the  k2  terms  are  to  retain  orthogon¬ 
ality.  The  kj  and  k2  are  arbitrary  gain  factors  which  will  be  discussed  later. 

It  appears  that  the  above  method  is  a  slight  improvement  on  the  version  given  by 
Corbett  in  that  the  orthogonality  correction  is  added  to  only  one  vector  rather  than 
to  both.  The  entire  drift  about  the  z  axis,  then,  is  determined  by  the  drift  of  i 
alone.  If  the  correction  were  fed  into  both  vectors,  the  total  drift  would  be  higher 
in  cases  where  the  two  vectors  tend  to  drift  in  the  same  direction. 

The  above  equatior.  ?  are  idealizations.  The  computer  which  is  used  to  instru¬ 
ment  them  will  be  solvii.  approximations  of  these  equations.  The  differences  give 
rise  to  errors  in  the  solutions  which  will  now  be  considered.  For  convenience, 
all  errors  are  divided  into  two  categories,  static  and  dynamic.  The  dynamic  errors 
are  associated  with  the  fact  that  the  actual  equations  the  computer  is  solving  are  of 
higher  order  than  the  ideal,  and  the  static  errors  arise  from  errors  in  resistance, 
capacitance,  pot  settings,  and  the  like.  These  two  types  of  errors  must  be  treated 
by  different  means.  The  dynamic  errors  will  be  considered  first. 

It  is  assumed  that  the  only  dynamic  effect  of  importance  is  the  bandpass  of 
the  multiplier.  Amplifiers  and  integrators  should  be  at  least  one  order  of  magnitude 
better  than  multipliers  in  this  respect,  so  the  assumption  appears  sound.  A  simpli¬ 
fied  analysis  will  illustrate  the  important  issues.  Consider  the  equations  for  com¬ 
ponents  of  the  vector  i,  with  correction  terms  deleted. 


a21R  "  a31^» 


a3 1 P 


lll 


R. 


(83) 


a31  =  allQ  ’  a2 1 P* 

Now  consider  a  special  -?vc.  -■  ere  P  =  R  =  0  and  Q  is  a  constant.  These 
vquatiore 
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all  =  a31°' 


a3i  =  anQ* 


a21  =  °* 


(84) 


If  a£j  is  initially  zero,  it  will  remain  so,  and  this  equation  may  be  deleted  from 
the  set.  If  Q  is  constant,  the  remaining  equations  are  linear.  Taking  the  Laplac 
transform  of  these  equations,  together  with  the  initial  conditions  a^j  =1,  a^j  ~  0, 
and  the  result  is 

(85) 

If  it  is  assumed  that  the  transfer  function  of  the  multipliers  is  G(S),  then  these 
equations  would  be 


Sall  =  1  "  a31°» 


Sa31  =  allQ 


Sa^j  =  1  -  a^QG(S), 

Sa^  —  anQG(S). 

It  is  possible  to  solve  these  equations  for  and  a^j. 


(86) 


11  "  S2  +  Q2G2(s)  ' 


QG(S) 

31  S2  +  QZGZ(S) 


(87) 


Before  proceeding  further,  it  is  necessary  to  make  some  assumption  concerning 
G(S).  It  should  be  clear  that  for  any  reasonable  result,  this  transfer  function 
should  be  only  slightly  different  from  unity.  It  must  equal  unity  at  very  low 
frequencies,  (static  errors  assumed  zero).  A  reasonable  assumption  is  that  its 
linear,  i.  e. ,  in  its  power  series  representation  the  non-linear  terms  are  negli¬ 
gible  in  comparison  to  the  linear  term.  We  have  assumed  as  reasonable 

G(S)  .It  TjS.  (88) 
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Substituting  this  into  Equations  (87)  gives 


(89) 


(90) 


f  = 

Several  conclusions  may  be  drawn  from  this.  If  r  is  positive  (corresponding 
to  a  lead  in  the  transfer  function)  then  the  length  will  decrease.  If  t  is  negative 
(corresponding  to  lag),  then  the  length  will  diverge.  The  term  t  Q  sin  2Qt  repre¬ 
sents  an  oscillatory  error  of  p  'ak  magnitude  tQ.  It  will  be  shown  later  that  the 
amplitude  correction  term  must  be  kept  as  small  as  possible  to  avoid  angular 
drift.  Therefore,  it  will  probably  not  be  possible  to  get  the. correction  gain  high 
enough  to  cut  down  this  oscillatory  error  term.  The  only  way  to  keep  it  small 
will  he  to  keep  tQ  small.  Let  us  assume,  for  instance,  that  a  computing  accur¬ 
acy  of  0.  1  per  cent  is  desired.  This  means  that  tQ  <  10"^.  It  may  also  be  seen 
that  tQ  is  simply  the  phase  lag  (in  radians)  at  the  frequency  Q.  For  the  example 
given,  the  phase  lag  at  the  frequency  of  oscillation  should  be  less  than  one  milli- 
radian,  or  about  0.  06  degrees.  The  correction  circuit  will  be  able  to  take  care 
of  the  long  term  exponential  increase  quite  well,  though  if  the  bandpass  requirement 
stated  above  is  met,  this  source  of  growth  of  the  vector  will  be  negligible  compared 
with  ■'  s-'  due  to  «tatir  err.-  .vhich  will  be  considered  next. 


•iQ!I 


t  Q  sin  2  Qt). 
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As  a  prelude  to  studying  the  effect  of  static  errors,  it  will  be  useful  to 
consider  a  general  symbolic  diagram  of  the  circuit  required  to  solve,  the 
equations.  This  is  because  the  problem  of  errors  is  inseparable  from  that  of 
scaling,  so  some  sort  of  scaling  must  be  assumed.  This  diagram  i.s  given  in 
Figure  1.  It  may  be  seen  that  there  are  36  multiplications,  six  of  which  are 
independent.  There  are  six  integrators,  six  summing  amplifiers,  and  in  any 
practical  circuit,  there  would  have  to  be  a  considerable  number  of  inverting 
and  isolation  amplifiers  as  well.  As  the  nature  and  number  of  these  depends 
on  characteristics  of  the  multipliers  being  used,  they  are  deleted  in  this 
figure.  A  practical  circuit  for  accomplishing  this  transformation  will  be 
considered  later.  It  is  convenient  to  introduce  a  special  notation  for  errors 

in  these  multipliers.  The  voltage  error  in  the  multiplier  which  multiplies  100a 

R  *  ^ 

and  100  ^ —  is  designated  ft  should  be  emphasized  that  this  is  the  actual 

m 

error  in  volts,  not  a  ratio  or  a  percentage. 

Now  assume  that  the  multipliers  have  infinite  bandpass,  but  do  have  static 
errors.  It  may  be  seen  from  the  circuit  of  Figure  1  that  the  equations  for  compo¬ 
nents  of  the  vector  i  are  actually 


kii  =  a21R-a3iQ+Kiaii(lX^i)+Rm^-Qm^iQ+^L^f 

i 

k21  =a3iP-auRtKlazl(iri)mm  («> 

'*31  ■  anQ  4ir-  -T5J- 

These  are  the  equations  which  the  machine  will  be  solving.  It  is  of  interest  to 
investigate  certain  properties  of  the  solutions  to  these  equations,  and  to  see  how 
they  compare  with  the  ideal, 

— ► 

First  consider  the  effect  on  the  length  of  the  unit  vector  i  .  Assume  that  Kj=  0 
so  that  no  corrections  are  neing  introduced.  The  equations  become 


all  '  a21R  '  a31°  +  Rr 


100  ~  °m  100 


a, .P-a,  ,R+  P  -TfTfT  '  R  “TTiTr  * 
31  11  m  100  in  100' 


allQ  "  a21P  +  Qm  100 


21 P 
100  * 
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Figure 


and,  difl  ci'cm  iatin  g  this 


* 


c* 


2  7  ?  Z 

The  length  of  i  is  given  by  ~  a  ^  ^  “  t  a^  j  "  +  a^j 
expression  with  respect  to  lime, 


Zii 


“al lall 


Za2ha21 


+  2. 


3 1 a  3 1 


The  length  is  close  to  unity  for  all  cases  of  interest,  so  that,  to  the  first  order. 


t 


allall 


a2la21 


a31a31‘ 


(94) 


Substituting  Equations  (93)  into  this  expression  gives 

;  E2IR  E31Q  ,  „  E31P  p 

11  m  100  al  i  m  100  a21  m  100  ■  a21  m  100 


hr 


(95) 


n  E11Q  p  E2 1 P 

a31Um  100  ‘  a31  m  100  * 


It  is  now  necessary  to  consider  the  nature  of  these  errors.  In  the  first  place, 
the  error  is  viewed  as  a  random  variable.  For  a  given  multiplier,  the  error  which 
exists  is  some  definite  function  of  the  two  inputs.  This  function  is  more  or  less 
repeatable,  at  least  over  a  short  time  period,  so  in  this  sense  it  is  not  a  random 
variable.  However,  when  one  considers  ciffeienl  multipliers,  the  error  existing 
at  certain  valu?s  of  the  inputs  is  now  a  function  of  which  multiplier  is  used  and  is 
thus  a  random  variable.  These  variables  are  considered  independent  because 
there  is  no  reason  to  expect  that  the  error  in  one  multiplier  will  influence  the 
error  in  another. 

Two  different  types  of  error  are  considered  in  the  subsequent  analysis.  In  the 
first  type,  called  "uncorrelated",  the  error  is  assumed  to  be  Gaussian  and  independ¬ 
ent  of  both  inputs.  In  other  words,  the  error  existing  at  any  given  value  of  the  inputst 

is  assumed  to  be  drawn  from  a  normal  distribution  of  zero  mean  and  variance  which  i 

« 

is  independent  of  either  input  and  of  the  particular  multiplier  being  used.  This  erroi; 
distribtuion  is  intended  to  be  consistent  with  electronic  multipliers,  though  no  propel1 
statistical  data  are  available  on  them.  Usually  electronic  multipliers  are  adjusted 

t 

so  that  the  error  when  either  of  the  inputs  are  zero  is  somewhat  smaller  than  other-, 
wise,  and  there  are  no  data  to  support  using  a  normal  distribution.  Until  such  data 
are  available  the  above  hyp-  theses  are  as  good  as  any  and  more  convenient  than  mos 

i 
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The  second  type  of  error,  called  "correlated"  differs  from  the  first  in  that 
the  variance  of  the  errors  is  assumed  proportional  to  one  of  the  inputs  and 
independent  of  the  other.  This  is  intended  to  represent  servomultipliers,  and 
the  experimental  justification  here  is  somewhat  better.  Since  this  is  the  type  of 
multiplier  used  in  the  simulation  described  later,  the  error  distribution  was 
measured.  Results  are  given  in  Appendix  C,  and  the  normality  of  the  distribution 
is  reasonably  well  verified.  Proportionality  of  the  variance  to  the  voltage  across 
the  potentiometer  was  not  checked  but  it  is  an  inevitable  consequence  of  the  nature 
of  servomultipliers. 

The  drift  rate  i  of  Equation  (95)  thus  becomes  a  random  variable.  It  is  the 
sum  of  six  independent  Gaussian  variables,  so  its  variance  is  the  sum  of  the 
variances  of  the  individual  terms.  Thus  the  standard  deviation  or  square  root  of 
the  variance  is  given  by 


*  =  100  fan2  Rm2  + 


a  2  P  2 
21  m 


<r.  is  the  standard  deviation  of  f  and  <r^  is  the  standard  deviation  of  multiplier 
errors.  It  is  not  possible  to  evaluate  this  expression  without  knowledge  of  the 
particular  Pm,  Qm  and  involved,  but  it  is  possible  to  determine  a  conveni¬ 
ent  upper  bound.  The  total  maximum  angular  rate  Wm  is  given  by 


W  = 
m 


'Q  2  +  P  2 
m  m 


Making  use  of  this  fact  together  with  the  normality  relation  aj  f  +  a22  +  a, 


it  may  be  seen  that 


In  the  foregoing,  errors  of  the  first,  or  uncorrelated  type  have  been  assumed. 

The  drift  rate  tends  to  be  proportional  to  the  full-scale  velocity,  regardless  of 
what  "velocity- actually  exists.-  -  -  •  ...  -  ...  ..  ...  . 

This  is  not  true  with  the  second  or  correlated  type  of  error.  Observe  that 
each  of  the  multipliers  is  used  to  multiply  a  direction  cosine  by  one  of  the  angular 
velocity  components,  P,  Q  or  R.  Consider,  for  example  the  multiplier  which 
generates  the  first  error  term  in  equation  (95);  the  one  multiplying  a^j  by  R.  . 
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If  the  direction  cosine  is  put  onto  the  shaft  of  the  scrvomult  ipli  e  r,  and  R 
is  put  across  the  potentiometer,  then  the  error  standard  deviation  will  be 
proportional  to  R  as  postulated  earlier.  This  fact  may  be  stated  as 


o 

t 


where  cr 

ch 


m 

is  the  standard  deviation  of  errors  at  full-scale,  i.  e.  ,  when  R  =  R 

rn 

other  multiplications  are  treated  similarly,  then  it  may  be  seen  that 


If  the 


cr 

I 


:  w  w 


(98) 


This  result  is  somewhat  more  favorable  than  that  for  the  uncorrelated  case. 

The  drift  rate  tends  to  be  proportional  to  the  existing  rate  of  rotation  rather  than 
the  full-scale  rate.  There  would  be.  little  difference  between  the  two  if  the  rate 
were  near  the  ma.\  mum  most  of  the  time,  but  this  is  not  usually  the  case.  In 
fact  in  many  simui,  -ions,  the  rotation  rate  is  small  most  of  the  time,  the  maximum 
rales  being  required  only  once  or  a  few  times. 

It  is  now  of  interest  to  investigate  the  effect  of  the  correction  scheme  on  this 
length  drift.  V'/o  consider  first  the  case  of  uncorrelaied  errors.  If  the  complete 
Equations  (92)  arc  substituted  in  Equation  (94).  the  result  is 


1  =  K.(l 


R  !A_ir  €jl£L  .  -  p  13J.P 

■11  tn  100  1 1  im  100  +  a21  m  100 


a  -  R  ~~  +  a  .  O  -yr—  -  a  P  }  (99) 

21  ni  1O0  31  m  100  31  in  100  J  v  ' 


K 


1  00C 


I 


j  1  J  1c  Z1  2)c  31  31c 


Consider  the  first  term  in  this  equation.  Observe  that  i*  i  is  simply  f  .  If  it 
is  as  s  timed  that  i  -  l  .0  +  Af ,  then  f  ^  -  1 .0  +  2Af,  so  that  the  first  term  becomes 
••  K.2AE  The  second  term  is  the  same  as  the  length  drift  rate  of  Equation  (95). 
The  last  term  is  the  drift  rate  due  to  the  correction  mechanism  itself.  For  the 
present,  it  is  assumed  that  is  large  enough  that  this,  last  tern:  is  negligible 

with  respect  to  the  first  two.  We  will  return  to  this  point  later.  vVith  this 
assumption,  then.  Equation  (99)  becomes 


„  c2 1 R  ^  *310 

-  .  H  -tPT'*  iQ  "T7V7T  ' 

1  1  m  1 00  11ml  00 


P  !ilp  Ills. 

21  m  100  "  21m  100 


(100) 


are 


C  >  -.lbiring  Equations  (  I  02)  anil  (103),  the  limits  on  K. 


\V  <y  W  C , 

_JIIJ - <  K,  £  — 

10o/2  <t  .  '  1  1 

At 


There  are  several  differences  for  correlated  errors.  When  P  =  Q  =  K  =  0, 
then  all  errors  will  be  zero,  provided  the  velocity  components  have  been  put 
across  the  multiplying  potentiometer s  as  suggested  above.  Under  these  circum¬ 
stances,  the  drift  is  determined  by  integrator  drift,  and  is  at  least  an  order  of 
magnitude  lower  than  the  drifts  normally  arising  from  multiplier  errors.  If, 
however  the  length  tolerance  is  to  be  met  under  the  largest  allowable  rates  of 
rotation,  then  the  lower  bound  on  is  the  same  as  that  of  Equation  (102). 

The  upper  b.  mid  will  disappear  if  the  error  signal  is  put  across  the  potentio¬ 
meter  in  the  (.  multiplications.  In  this  case,  the  error  arising  in  the  multi¬ 
plier  would  be  o  the  order  of  rather  than  of  the  order  of  «  .  The  gain  would 
not  have  an  upper  bound  then,  except  lor  the  fact  that  integrators  tend  to  drift 
s oil* g wlidt,  faster  tlic  1  * i c r  chains  prccceclii"*^  »hcin. 

Next  the  angular  drift  rate  will  be  considered.  Even  if  the  coordinate  system 
is  kept  orthogonal,  there  is  still  the  tendency  lo  drift  in  orientation,  which  comes 
principally  from  static  errors  in  the  multipliers.  The  components  P,  Q  and  R 
express  the  angular  velocity  in  the  moving  system.  If  P  ,  Qq  and  are  the 

components  in  the  fixed  system,  it  m;  _•  be  shown  that 


Ro  “  al 1 3 1  2  +  3 2 1 n 2 2  f  a31a3Z  J 


Ro  ~  3 1  2 3 1  1  f  a22a21  +  a32a31  ' 


P  +anar  +  a.  -.a  7  , 

o  1  Z  1  3  Z  Z  Z 3  3  Z  3  3 


(101) 


P  =  a.,  a.,  +  a  ,  a  +  a ,  .  n  , 

o  1  3  1  2  2  3  22  3  o  32 


Qo  "  al 3 3 1 1  +  3 2 3 3 2 1  +  3 3 3a 3 1  ’ 


Q  =  a .  .  a .  .  +  a  a  +  a  ,  .  a 
o  1 1  1  J  21  23  31  33 


'.V  : 


'<  *<- 17 
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This  is  a  first  order  differential  equation  for  A f,  the  right-hand  side  or  forcing 
term  being  a  combination  of  errors.  As  the  errors  change  during  the  course  of 
a  run,  then  the  length  error  A /  will  change  also,  but  for  simplicity  consider  a 
static  case.  If  the  coordinate  system  is  not  rotating,  then  all  direction  cosines 
are  fixed,  and  the  right  side  of  Equation  (  100)  becomes  a  constant.  Let  us  call 
it  u.  This  constant  will  have  different  values  depending  on  the  set  of  multipliers 
used,  and  in  fact  the  standard  deviation  of  the  values  it  can  assume  is  given  by 
Equation  (97).  In  the  steady-state,  Ai  must  be  zero,  and  from  Equation  (100) 


2Kj  Ai  =  u 


The  deviation  in  Ai  will  then  be  given  by 


°Ai  =  “ZE7  ®u  = 


KjIOO/T  m  * 


(101) 


<r4  is  fixed  by  the  nature  of  the  multipliers,  tfmby  the  scaling  requirements.  Kj 
may  be  chosen,  however,  to  make  9 .  j  as  small  as  desired.  Conversely,  the 
required  value  of  Kj  is 


K1  > 


ioo  ST  i 


(102) 


It  was  mentioned  that  the  third  term  of  Equation  (99)  should  be  negligible  with 
respect  to  the  first  two.  The  reason  for  this  is  not  clear  from  what  has  been 
presented  thus  far,  but  it  will  be  seen  later  when  angular  drift  is  considered. 

Angular  drift  is  determined  by  a  similar  equation,  and  it  seems  unreasonable  to 
allow  the  length  correction  circuit  to  contribute  to  the  drift  in  angle  when  it  can  be 
avoided,  as  will  now  be  shown,  by  proper  distribution  of  gains.  The  third  term  of 
Equation  (99)  will  have  the  standard  deviation  j  <r  ,  while  the  standard  deviation 

of  the  second  is  If  it  is  required  that  the^third  term  be  less  than  one 

quarter  of  the  second,  then 


100C, 


(103) 


K1  /Wm 
C1  4 

which  determines  C^,  once  is  chosen.  This  establishes  the  upper  limit* on  Kj, 
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Now  substituting  Equations  (92)  into  the  third  of  these  gives 

O  o  ^  r>  .  f  „  e21R  e31Q 

-R  =  -a.P-a  „Q.a,)Rt{al)R  — .  a  yr-  -  a.  ,Q  — ry,n 
o  13  23  33  1  u  m  100  12  m  100 


4  a  P  R  IllJl  n  !i'_Q  ‘  p  l^JP- 

22  m  100  "  a22  m  100  '  a32Um  100"  '  a32  m  100  ■ 


Kj 

lOOC^  a12  *110  +  a22  *2 1  c  +  a32  *31ci  *  (105) 

The  first  three  terms  represent  the  transformation  of  the  velocity  vector  into  the 
fixed  axis  system.  In  other  words,  they  represent  the  "correct"  value  of  Rq.  The 
remaining  terms  represent  the  error  in  RqI  or  the  component  of  the  drift  velocity 
vector  along  the  unprimed  Z  axis.  Combination  of  the  errors  is  similar  to  that 
for  drift  in  length.  If  it  is  required  that  the  last  term  of  Equation  (105)  be  negligible 
with  respect  to  the  others,  then  the  standard  deviation  of  drift  rate  is  given  by 


W 


m 


"c  U 
100 


V/z 

100 


(106) 


for  the  uncorrelated  and  correlated  error  cases  respectively.  Thus  it  appears  that 
the  drift  rate  will  be  a  constant  fraction  of  full  scale,  for  uncorrelated  errors,  and 
will  be  a  constant  fraction  of  the  angular  velocity  for  correlated  errors.  The 
conditions  on  Kj  and  Cj  are  the  same  as  developed  earlier,  namely 


cr 

t _ 

100  fl  ^ 


(107) 


The  drifts  in  P  and  Q  are  substantially  the  same  provided  gains  are  chosen 
such  as  to  make  the  drift  contribution  of  the  correction  circuits  negligible.  The 
case  of  and  is  substantially  the  same  as  for  and  Cp  and  they  should  be 
chosen  in  the  same  way. 

Analyses  of  drifts  in  P  and  Q  are  done  is  much  the  same  way  and  lead  to 
similar  results. 
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B.  Quaternion  Error  Analysis 

The  quaternion  simulation  may  be  handled  in  much  the  ssmfi  way.  The 
quaternion  components  are  not  all  independent,  and  v,c  may  make  use  of  the 
relationship, 

ej2  +  ez2  +  e32  +  e44  =  1  ,  (32) 

in  the  same  way  that  +  a21^  +  a31^"  =  *  was  used  to  maintain  the  length 

of  the  unit  vector.  The  equations  to  be  solved  are 


2ej  =  -  e^P  -  e^Q  -  e,R  +  Kij>e.j  , 


2e^  =  -  e^P  +  e^Q  +  e^R  +  Ki^e,  , 


2e 


3 


e^P  *  el^  "  eqR  +  K4>e3  » 


(108) 


2e^  =  e.  P  -  e^Q  +  e^R  +  Kbe^  , 


<t>  = 


Hie  bandpass  rrquiremm  oi  this  metnod  is  only  one  nail  as  severe  as  that  for 
direction  cosines.  This  may  be  seen  in  several  ways.  Consider  the  equations 
above  with  Q  constant  and  P  =  R  =  K  =  0. 


c 


e 


e 


c ,  Q  , 
2 


2  =  + 


3  =  +  el-Q-> 


-  e,  _Q 
4  2* 


(109) 


For  initial 
conditions. 


conditions,  wc  assume  e 
e2  ar'4  '-i  '-'/ill  remain 


l 


=  1, 

zero 


e2  ~  e3  “  e4  =  0.  Under  these 
and  the  equations  become 


=  -  « ,  Q.  , 


ei  jQ_ 
*  2 


(no) 
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These  are  identical  with  thoae  treated  in  the  bandpass  study  o£  the  proceeding 
Section  except  that  the  frequency  ie  reduced  by  one  half.  The  tame  conclusion 
can  be  seen  from  the  definition  of  «j,  From  liquation  (33),  e^  =  cos  p/2,  so 
while  u  is  completing  a  full  rotation,  p/2  only  moves  through  180°.  This  means 
that  the  servos  which  are  driven  by  the  e's  only  move  with  half  the  speed  of  thoae 
driven  by  the  a's,  for  a  given  rate  of  rotation  of  the  coordinate  system.  Conse¬ 
quently,  for  a  given  accuracy,  twice  as  much  phase  shift  may  be  allowed. 

In  order  to  analyze  the  effect  ol  static  multiplier  errors,  it  is  again  necessary 
to  postulate  some  particular  scaling.  The  simplified  diagram  is  shown  in  Figure  2. 
The  equipment  necessary  for  determining  the  direction  cosines  has  been  included, 
as  they  will  always  be  needed.  The  basic  quaternion  component  computation  requires 
20  multipliers  and  four  integrators.  Conversion  to  direction  cosines  requires 
another  six  multipliers,  for  a  total  of  26  multipliers  and  4  integrators,  against  36 
multipliers  and  6  integrators  for  direction  cosines.  Notation  for  individual  multi¬ 
plier  errors  is  similar  to  that  applied  earlier.  The  voltage  error  in  the  multiplier 
which  multiplies  R  and  e,,  for  instance,  is  called  «^,  and  soon.  From  this 
figure,  it  may  be  seen  that  the  equations  being  solved  are: 

.  o  F 

m  iuu 


4e 


1 


_  •  ft  -  O  X  VAa 

c4  "  °3W  '  C24V  '  *v^l  ~ 


'4P  100  '  ~m  iuu 


I  M.*,  ~  _  >  AA 


2e 


2  =  -  e^P  +  e^Q  +  e^R  +K<|»e^ 


*3P 


m  100'+  Qm  100 


1±Q-+ 

i  no  ’ 


*1R  K  *2c 
ml00+  c  100 


(HI) 


.  _  P 

2e3  =  +  F  e j  Q  -  e^R  +  K<j>e^  +  Pm'Joo 


2e.  =  +  e,P  -  e,Q  +  e?R  +  K*e4  +  PmT 


1P 


00 


..  n  R  Sr  ,  K  *3c 

Sn  100  "  m  100  c  100 


x  +  ti  I?*  .  K  *4c 

V  100  m  100  c  100 


Consider  first  the  effect  of  errors  on  the  length  of  the  quaternion.  The  length 

is  given  by  4  +  e^  +  e^  .  Differentiating  with  respect  to  time  give* 

•  •  •  •  | 

if  *  •]*!  •+■  ®2e2  +  e3°j  +  -Since  the  length  is  not  to  be  allowed  to  vary  much 

from  unity,  the  drift  rate  may  be  approximated  by  f  =  ej®i  +  +  e3®3  +  e4®4* 

Substituting  Equations  (ill)  into  this  relation  and  performing  the  appropriate 
reduction*  gives: 
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(112) 


If  /  =  1  +  A/,  then  the  first  term  is  simply  -  2KAI.  Thus  Equation  (112)  becomes 
A1  +  KAf  =  Jqq  {  -el  Pm«4p  -  eiQm  *3Q  *  eiRm  «2R  '  *2Pm  *3P 


+e2  Qm  *4Q  +  e2Rm  *1R  +  e3Pm  *2P  +  e3°m  CIQ 


e-R  < . D  +  e.P  -  e.Q  t.  +  e.R^  *,„} 
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+  2000^  el  ‘lc  +  *2  *2c  +  *3  *3c  +  ®4  *4c  ^ 


(113) 


which  is  analogous  with  Equation  (100). 

As  in  the  proceeding  section,,  we  assume  the  second  term  on  the  right  of 
Equation  (113)  to  be  negligible  with  respect  to  the  first.  The  variance  of  the 


first  term  is 


Thus,  for  the  steady-state. 
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where  <r.  .  is  the  standard  deviation  in  the  length  erroi  which  will  be  allowed. 

A *  iC  . 

The  standard  deviation  of  the  second  terms  of  Equation  (113)  is  Yqqq  tr<«  The 

requirement  that  the  second  term  be  no  more  than  one  quarter  the  first  gives 

K  C 
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Combining  the  two  gives  almost  the  same  conditions  on  K  as  were  obtained 
for  xbe  direction  cosine  method. 


Wr 

200  <r. , 


.  ^  W 

<K  C  . 


(115) 


The  foregoing  is  for  uncorrelated  errors.  The  changes  in  the  correlated  error 
case  would  be  the  same  as  in  the  direction  cosine  analysis:  the  upper  bound  on 
K  disappears  and  the  lower  bound  is  the  same  as  in  the  uncorrelated  case. 

Next  the  effect  of  multiplier  errors  on  the  angular  drift  must  be  evaluated. 
Equation  (79)  gives  the  rotation  rata  components  P,  Q,  and  R  in  terms  of  the 
rates  of  change  of  the  quaternion  components.  The  apparent  value  of  P,  which  we 
we  call  P^  is,  the  P^  =  2(-e^e^  -  e^e^  +  e^e^  *  ®1®4^*  ^  Equations  (111)  are 
substituted  into  this  expression,  the  result  is 

P.  *  <V 2  ♦  L  «4*>  P  *  TOT  WPm  «4P  ♦  °n,  *JQ  +  Rm  «2R> 

+  ®3^Pm  *3P  '  Qm  *4Q  "  Rm  ‘lR*  +  e2^Pm  *2P  +  Qm  *1Q  “  Rm  *4R*  ^ l6^ 

+  ®l*Pm  ‘lP  ‘  Qm  *2Q  +  Rm  ®3R*  }  +  rfe^"®4  ‘lc  ‘®3*2c  +®2  *3c  +  V4J  * 


The  first  term  is  equivalent  to  (  1  +  2Af)Pso  the  error  in  P,  (P  -  P)  is 

A 

AP  =  Pa-P  »  2Afp  +  Too{«4(pm«4p  +  Qm  *3Q  +  Rm  «2r5 


(117) 


+  ®3^Pm  *3P  "  Qm  *4Q  "  Rm  ‘lR*  +  ®2  ^Pm  *2P  +  Qm  *1Q  ‘  Rm  *4R* 

+  ®l^Pm  *1P  “  °m  *2Q  +  Rm  *3R*  }  +  To5c{’®4*lc"  ®3*2c  +  e2*3c  +  ®1  *4c^‘ 

2Af  is  related  to  the  choice  of  K.  From  Equation  (114)  it  may  be  seen  that  the 

standard  deviation  in  Af  is  given  by  g  2~o*OK  before,  we  require  the 

last  term  of  Equation  (117)  to  be  negligible.  It  remains  to  determine  the  variance 

of  the  second  term.  By  taking  the  sum  of  variances  of  the  individual  terms,  it  is 

possible  to  show  that  the  standard  deviation  of  the  second  term  is  W  <r  . 

iuu  m  1 
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Thus  the  standard  deviation  in  AP  is 


AP 


100 
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(118) 


This  is  comparable  to  the  value  obtained  for  direction  cosines  if  K  is  made 
reasonably  large. 

Exactly  analogous  results  are  obtained  for  errors  in  the  other  two  axes. 
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Figure 


SECTION  VII 


SIMULATOR  RESULTS 


The  next  step  in  investigation  of  quaternion  coordinate  conversion  was 

solution  of  the  equations  on  an  analog  computer.  The  direction  cosine  method 

was  similarly  investigated,  and  an  attempt  was  made  to  make  the  conditions  of 

the  two  simulations  as  nearly  alike  as  their  inherent  differences  would  permit. 

Both  were  done  on  REAC  Series  100  equipment.  ServornulLplier s  equipped  with 

potentiometers  of  0.05  percent  linearity  were  used  throughout.  The  multipliers 

were  not  specially  calibrated  for  these  simulations,  so  their  adjustment  was 

consistent  with  normal  practice  in  the  Analog  Computation  Branch,  Aeronautical 

Research  Laboratory.  Of  course,  correct  pot  loading  was  used  in  all  cases.  In 

both  simulations,  a  maximum  rotational  velocity  of  0.  5  radians/sec  was  employed 

though  scaling  was  such  that  P  =  Q  =  R  =1.0  and  both  directions  cosines 

°  m  m  m 

and  the  quaternion  components  were  scaled  so  that  50  volts  (out  of  100  full  scale) 
represented  the  extreme  possible  excursion  of  the  variable.  This  was  done  to  get 
away  from  possible  end  effects  on  s ervomuitiplier  potentiometers. 

Adequate  checking  of  a  coordinate  conversion  is  a  problem  in  itself,  and  while 
it  is  not  claimed  that  the  method  adopted  lien;  meets  til  requirements,  it  seems 
sufficient,  and  no  better  method  presented  Rself.  With  the  types  of  coordinate 
conversion  considered  here,  there  are  two  things  principally  to  be  checked:  the 
action  of  the  "orthogonalization"  or  correction  mechanisms,  and  the  rotation  drift 
rate.  The  first  may  be  checked  simply  by  monitoring  the  error  quantities  which 
are  used  for  correction.  The  drift  rate  is  not  so  easily  checked.  In  most  cases 
the  drift  rate  will  be  very  small  compared  with  the  rate  of  rotation.  An  exception 
to  this  is  the  case  when  an  angular  velocity  of  zero  is  desired.  Any  shift  which 
takes  place  under  these  conditions  is  readily  detected.  To  check  the  drift  while 
rot  iting,  the  following  procedure  was  adopted:  a  single  input  of  P  =  0.  500  rad/sec 
(O  and  R  zero)  was  applied  to  the  equations  for  a  period  of  approximately  125.66 
seconds.  This  is  enough  time  for  ten  complete  revolutions  at  this  frequency.  The 
transformation  matrix  should  be  the  same  at  the  end  of  this  period  as  it  was  at  the 
beginning,  except  for  drift.  At  the  start  of  each  run,  the  transformation  was  the 
identity  transformation,  so  the  matrix  existing  at  the  end  of  the  run  has  a  simple 
interpretation.  Because  of  the  difficulty  of  controlling  the  length  of  run  with 
sufficient  accuracy,  several  runs  were  made  using  the  same  conditions,  and  the 
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time  of  run  was  recorded.  It  was  then  possible  to  plot  the  rotation  angle 

1;  Junction  of  run  n-r.e,  and  by  interpolation,  to  determine  the  drift  angle 
»  m.'.u  at  exactly  1  d 5 ..  06  sei-onus.  The  same  process  was  repeated  for 

p  _  p  -  c,  Q  =  0.  5  and  P  =  Q  =  0,  R  -  0.  5,  thus  giving  rotation  about  each  of 

the  three  axes  singly. 

A.  The  Quaternion  Method 

It  was  mentioned  earlier  that  the  problem  was  scaled  for  50  volts  maximum 
on  the  e’s  rather  than  100  volts.  For  this  case.  Equation  (118)  becomes 


tt  W 
50  m 


fZ  K2 


(119) 


The  value  of  P  was  0.  5,  K  was  2.0,  and  servcmultipliers  were  used,  so  for 
correlated  errors,  this  expression  becomes 


“ap  ~  100  O*03)* 


(120) 


It  is  shown  in  Appendix  C  that  a'  is  0.05  volts.  Therefore,  the  standard 

-4  * 

deviation  in  P  is  5.  1  5  x  1  0  "*  raJiar.s/soccnd,  In  123  seconds,  this  would 
amount  to  0.  0644  radians  or  about  3.7  degrees.  The  standard  deviation  of  the 
drift  angle  after  125  seconds,  then  should  be  about  3.  7  degrees.  The  drift  angle 
was  determined  only  three  times,  once  each  for  P,  Q,  and  R  inputs.  Three 
results  constitute  an  insufficient  number  of  sample*  for  statistical  significance. 

In  order  to  get  the  number  of  results  required,  it  would  be  necessary  to  do  the 
entire  setup  many  times  using  different  pots  in  different  permutations.  It  was 
not  fjclt  that  the  improved  confidence  in  the  error  analysis  would  justify  the 
irmneme  labor  of  this  procedure.  The  results  of  the  three  determinations  which 
v.ero  r.’dc  are  not  inconsistent  with  the  theoretical  errors  found. 

At  the  end  of  each  run,  the  transformation  matrix  existing  is  very  ncar’y 
the  identity  transformation.  In  order  to  interpret  this  final  matrix,  it  is  conven¬ 
ient  to  make  use  of  Equation  (l  6  1)  of  Appendix  A, 


(A)  = 


0  1  0 
0  P  1 


cos  P 


Ap  cos  y 


-Ap  cos  a 


-A p  cos  p 
Ap  cos  a 
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drift  angle,  and  a,  p,  and  y  are  the  angles  between  the  drift 


axis  ami  the  x,  y,  and  /,  ax o s  respectively.  Thus,  after  the  run,  ajj»  & 2.1 
and  should  be  um'y.  All  the  other  direction  cosines  should  be  small,  but 

may  differ  front  zero.  If  they  are  small,  the  following  relationships  should 
hold: 

ald  =  "a21*  a13  =  ■a31*  aZ3  =  “a32* 

The  procedure  was  as  follows.  At  the  start  of  each  run,  initial  conditions 
of  the  four  integrators  were  set  to  correspond  to  the  identity  transformation, 
that  is  Cj  =  1,  e^  -  0,  =0,  =  0.  Then  the  computer  was  put  into  the 

"Operate"  condition.  This  switching  also  started  a  timer  driven  by  a  synchro¬ 
nous  motor  from  the  60  cycle  line  voltage.  An  integrator  and  biased  relay  were 
used  to  terminate  the  run,  and  this  termination  stopped  the  timer.  The  time  of 
the  run  could  then  be  read  directly.  At  termination,  the  computer  was  put  into 
"Hold"  and  voltages  corresponding  to  the  nine  direction  cosines  were  read  to 
the  nearest.  10  millivolts  with  a  digital  voltmeter.  The  results  obtained  are 
included  in  Table  I.  In  this  table  also  are  included  the  drift  angle,  and  direction 
cosines  of  the  drift  axes  for  each  run  made. 
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drift  runs,  quaternion  method 


The  first  five  runs,  those  made  with  P  =  0.  5,  Q  =  R  =  0,  may  be  used 
to  illustrate  several  points  of  interest.  In  all  eases,  the  drift  angle  will 
be  small,  corresponding  to  an  infinitesm.il  rotation,  .Since  an  infini'  .-.mal 
rotation  may  be  treated  as  a  vector,  we  may  take  components  of  this  vector 
along  the  three  axes  of  the  moving  system;  these  components  are  simply 
the  ^23*  a3J  and  aj2  matrix.  These  are  plotted  as  functions  of 

time  in  Figure  3.  For  this  case,  perfect  performance  would  dic  tate  that  the 
y  and  z  components  remain  zero,  and  the  x  component  increases  at  a  rate 
of  2 8,65  dcg/sec  (C.  5  rad/sec),  passing  through  zero  at  125.66  seconds. 

This  latter  is  shown  as  the  heavy  solid  curve  of  Figure  3.  It  may  readily  be 
seen  how  the  results  achieved  vary  from  this  simplification.  The  z  error 
is  about  0.2°,  the  y  error  is  about  0.6°  and  the  x  error  is  about  5,35°. 

This  illustrates  the  result  that  in  every  case  the  major  portion  of  the  drift  is 
in  the  direction  of  rotation.  This  is  to  be  expected  when  servomultipliers 
ar  used. 

Figure  4  shows  the  total  drift  angle  as  a  function  of  time  for  all  three 
cases.  Again  the  ideal  curve  is  the  heavy  line.  The  P,  Q  and  R  cases  show 
5.2,  7.0  and  3,4  degrees  respectively.  It  appears  also,  that  all  three  curves 
lie  significantly  below  the  ideal.  This  was  not  predicted  in  the  error  analysis. 
The  bias  appeared  in  the  direction  cosine  simulation  as  well,  and  a  more  de¬ 
tailed  consideration  will  be  given  to  it  in  the  next  part.  For  the  present,  it 
will  simply  be  said  that  it  was  traced  to  the  fact  that  the  gains  of  the  inverting 
ampiifiers  in  the  KEAC  are  consistently  slightly  less  than  the  indicated  value. 
Aside  from  the  bias,  it  seems  that  the  dispersion  agrees  relatively  well  with 
the  predicted  standard  deviation  of  3.  7  degrees.  This  does  not  mean  much, 
however,  in  the  presence  of  the  bias.  There  is  no  reason  to  expect  that  the 
bias  will  he  exactly  the  same  in  all  three  cases,  so  it  cannot  be  determined 
what  part  of  the  dispersion  is  due  to  multiplier  inaccuracy,  and  what  is  due  to 
amplifier  gain  variations. 

It  was  found  possible  ir.  the  direction  ccsine  method  to  reduce  the  bias 
markedly  by-  trimmi..g  the  amplifier  gains  to  exactly  the  desired  value.  -  Ir  is 
felt,  however,  that  the  results  should  be  presented  as  originally  obtained, 
however,  because  of  the  avowed  objective  of  showing  what  might  be  obtained 
in  a  practical  simulation  program.  If  U  is  desired,  the  drift  could  no  doubt  be 
reduced  by  detailed  calibration  and  adjustment  to  about  one  tenth  of  that  shown 
in  Figure  4, 
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TIME.  (SECONDS} 


Figure  5  shows  the  action  of  the  correction  circuit  as  it  eliminates  an 
5" :iial  error. 


Figure  5.  Correction  Circuit  Effect 

Figure  6  shows  a  Mme  history  of  part  of  one  of  the  runs.  It  shows  how  the 
correction  circuit  main  ains  the  length  during  a  run.  The  indicated  At  is  less 
than  0.  1  per  cent  most  oi  the  time.  This,  of  course,  is  merely  the  amplitude 
of  the  error  signal.  It  does  not  necessarily  mean  that  t  is  actually  being  held 
to  this  absolute  precision,  but  rather  that  it  is  being  held  to  the  value  of  1 
which  the  computer  shows  to  be  unity. 

The  drift  when  P  =  Q=  R  =  0isof  considerable  interest.  In  many  simu¬ 
lations  which  make  use  of  a  coordinate  conversion,  the  rate  of  rotation  is  small 
most  of  the  time,  reaching  the  peak  values  only  occasionally.  In  such  cases, 
the  tendency  to  drift  when  the  coordinate  system  should  be  standing  still  is  of 
prime  importance.  For  uncorrelated  errors,  there  is  no  reason  to  expect  that 
the  drift  will  be  much  smaller  in  this  condition  than  otherwise.  It  is  true  that 
the  errors  in  electronic  multipliers  are  somewhat  smaller  near  zero  than  else¬ 
where,  but  the  difference  is  not  dramatic.  This  is  particularly  true  when  one 
of  the  inputs  is  large  and  the  other  is  zero.  This  case  will  arise  since,  regard¬ 
less  of  the  orientation  of  the  coordinate  system,  some  of  the  e's  will  be  large. 
For  servomultiplicrs,  on  the  other  hand,  the  error  tends  to  zero  as  the  voltage 
across  the  multiplying  pot  goes  to  zero.  This  would  seem  to  show  that  the 
drift  would  be  exactly  zero  with  no  input  rate.  This  is  not  correct,  of  course, 
because  integrator  drift  is  still  present.  Rurts  18  through  27  of  Table  1  were 
made  under  these  conditions.  Rimr  '.3  through  20  were  n<*  -e  on  one  day, 
and  a  drift  rate  of  7.  ,c  •'  ^eg/sec  was  observed.  t'his  amounts  to  about 
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0.  0128  per  cent  of  full  scale  (1  rad/sec).  Runs  21  through  26  were  made 
about  a  week  after  the  preceeding  set,  after  another  problem  had  been  on 
the  computer  in  the  meantime.  From  this  second  set,  a  drift  rate  of  about 
2.  5  x  IQ'3  deg/sec  was  determined,  an  improvement  by  about  the  factor  3. 
Integrators  had  been' balanced  in  both  cases,  so  the  only  conclusion  possible 
is  that  the  zero -input  drift  is  somewhat  variable.  It  was  found  that  the  zero- 
input  drift  was  proportional  to  the  maximum  rate  for  which  the  computer  is 
scaled.  The  above  results  were  taken  for  1  rad/sec  full  scale,  but  some 
runs  were  made  with  5  rad/sec  full  scale,  and  the  drift  was  almost  exactly 
five  times  as  great. 

While  it  does  not  appear  utile  to  give  the  complete  computer  diagram, 
some  remarks  concerning  the  setup  are  in  order.  In  order  to  accomplish  the 
functions  indicated  in  Figure  2,  it  was  found  necessary  to  use  36  summing 
and  inverting  amplifiers,  4  integrators,  and  8  multiplying  serves,  each  with 
three  multiplying  pots.  One  summing  amplifier  was  used  ahead  of  each  inte¬ 
grator  to  do  the  summing.  now  appears  that  this  was  not  wise,  because 
while  the  integrator  gains  proved  quite  accurate,  the  summing  amplifier  gains 
were  less  so,  and  introduced  an  error  in  angular  rate.  The  simulation  was 
not  found  to  be  critical  or  sensitive  in  any  way,  except  that  the  zero -input 
drift  varies  somewhat  from  day  to  day.  In  all  cases,  however,  it  was  quite 
acceptable. 

B.  The  Direction  Cosine  Method 

The  direction  cosine  simulation  was  done  under  as  nearly  the  same  conditions 
as  pcosible.  The  same  rate  of  rotation  (0.  5  rad/sec)  and  the  same  full-scale 
rate  were  made  and  the  same  length  of  time  was  used.  It  was  found  necessary 
to  time  the  runs  with  a  Berkley  counter,  rather  than  the  syncrhonous  clock  used 
in  the  quaternion  case  because  of  the  failure  of  the  latter  unit  at  the  start  of  the 
cvame  runs.  Such  checks  as  could  be  made  showed  no  significant  difference 
in  the  timers.  To  do  the  operations  of  Figure  1,  31  inverting  and  summing 
amplifiers,  6  integrators  and  12  s ervomultipliers  were  required.  Table  2  shows 
tne  results  obtained  from  reading  the  direction  cosines  at  the  end  of  each  run. 
Again  the  rotation  angle  was  plotted  against  time  for  each  of  the  three  inputs. 
These  plots  are  given  in  Figure  7.  Ap’;n  it  appears  that  the  three  curves  show 
something  like  the  predict--;  di=*~  rsion,  but  a  large  bias  also.  The  bias  in  this 
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case  is  nearly  twice  that  found  for  the  quaternions.  The  cause  proved  some- 
what  difficult  to  locate,  and  it  is  constructive  to  consider  some  of  the  checks 
that  were  made  in  the  process  of  looking  for  it. 

This  bias  amounts  to  an  erroneous  rate  of  rotation.  The  coordinate  system 
is  rotating  a  trifle  too  slowly.  In  the  example  considered  here,  the  system  has 
completed  about  10  revolutions,  turning  through  a  total  angle  of  nearly  3600 
degrees.  At  the  end  of  this  time,  it  is  in  error  by  some  10  degrees  or  so.  This 
is  one  part  of  360,  surely  not  a  large  error,  yet  it  is  the  predominant  one,  being 
nearly  three  times  as  large  as  the  error  due  to  multiplier  static  errors.  This 
same  bias,  though  smaller  in  size  was  observed  in  the  quaternion  simulation. 
This  leads  to  a  strong  suspicion  that  it  is  due  to  some  characteristic  of  the 
computing  equipment  itself,  rather  than  some  outright  mistake  in  the  setup. 

To  be  sure,  the  setup  was  checked  most  carefully,  many  components  were 
interchanged  and  all  the  usual  checking  methods  were  applied  with  no  result 
other  than  to  reaffirm  that  the  equipment  was  correctly  wired  in  accordance 
both  with  the  circuit  diagram  and  the  equations.  Then  checking  of  the  counter 
characteristics  was  started.  It  was  thought  that  possibly  the  time  base  of  the 
computer  (i.  e.  the  time  constant  of  the  integrators)  was  not  exactly  the  same 
as  that  of  the  counter.  To  check  this  a  linear  0.  5  rad/sec  oscillator  of  50  volt 
amplitude  was  set  up  using  only  two  integrators,  two  hand -set  potentiometers 
and  one  inverting  amplifier.  This  oscillator  was  then  allowed  to  run  for  125.66 
seconds,  and  the  result  measured.  The  oscillator  agreed  precisely  with  the 
timer.  The  output  of  the  oscillator  integrator  which  had  a  zero  initial  condition 
(the  other  one  had  a  50  volt  initial  condition)  went  through  zero  within  10  milli¬ 
seconds  of  125.66,  and  the  time  error  we  are  looking  for  is  more  in  the  nature 
of  300  milliseconds.  This  is  a  surprisingly  good  check. 

Next  it  was  thought  that  it  might  be  a  phase  lag  effect.  This  was  ruled  out 
by  two  experiments.  First,  a  run  was  made  so  as  to  allow  only  five  complete 
rotations  of  the  system  rather  than  10.  It  was  found  that  the  bias  was  very 
nearly  one  half  that  of  Figure  7.  If  it  had  been  a  phase  lag  effect,  it  would  have 
"been  more  nearly  constant  with  time.  Furthermore,' the  phase  lag-of  a  m-ulti--  - 
plier  was  measured  directly  at  this  frequency,  and  was  found  to  be  about  0.  6 
milliradians.  The  phase  shift  required  to  explain  this  bias  would  be  more  in 
the  order  of  0. 1  radians. 

By  setting  up  a  separate  oscillator  using  two  integrators  and  two  servos, 
the  difficulty  was  finally  traced  to  the  fact  that  the  summing  and  inverting 
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amplifier  gains  are  consistently  low.  The  gains  were  carefully  adjusted  in  this 
c  ..  j. irate  oscillator,  and  most  of  the  bias  was  removed.  It  was  not  considered 
•  orthwhile  to  similarly  trim  all  the  summers  and  inverters  in  the  entire  simu¬ 
lation. 

After  all,  the  aim  of  this  program  was  to  show  typical  simulator  results. 
Still,  if  the  ultimate  accuracy  of  which  the  servos  are  capable  is  to  be  attained, 
something  of  this  nature  should  be  done.  It  should  be- mentioned  that  not  all  of 
the  amplifiers  would  have  to  be  trimmed.  Only  those  amplifiers  which  are 
between  the  P,  Q,  and  R  multiplications  and  integrator  inputs  are  critical. 

None  of  those  in  the  correction  loops  can  give  trouble.  This  is  true  of  both 
methods.  It  is  rather  surprising  to  find  that  servos,  generally  viewed  with 
suspicion  and  avoided  when  possible,  should  not  prove  to  be  the  major  source 
of  error  in  these  simulations.  Amplifiers,  rather,  have  proved  to  be  the 
limiting  factor. 

It  may  also  be  showi  from  the  data  of  Table  2,  that  the  zero -input  drift  of 
this  method,  at  least  on  t,  ?  day  the  runs  were  made,  was  about  3  x  10"^  deg/sec. 
This  is  comparable  with  the  better  set  of  results  obtained  with  quaternions. 
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SECTION  VIII 


SUMMARY  AND  CONCLUSIONS 

< 

Having  investigated  both  quaternion  and  direction  cosine  coordinate 
conversions,  it  is  now  in  order  to  compare  the  strengths  and  weaknesses 
of  the  two  methods.  For  comparison  purposes,  the  Euler  angle  method 
will  be  included,  though,  as  mentioned  earlier,  it  is  not  str:ctly  comparable 
with  the  other  two  in  capability.  Before  proceeding  with  this  comparison,  it 
is  well  to  outline  the  criteria.  It  seems  that  there  arc  two  main  areas  of 
comparison.  Listed  in  order  of  their  importance  they  are  (l)  minimum 
equipment  requirement  (both  qualitative  and  quantitative)  and  (2)  ease  of 
programming,  reliability  and  comprehensibility. 

As  far  as  qualitative  requirements  on  equipment  is  concerned,  it  has 
been  shown  that  both  direction  cosine  and  quaternion  methods  a:e  equally 
sensitive  to  multiplier  errors.  It  further  seems  clear  that  since  the  Euler 
angle  system  (as  proposed  by  Howe  and  Gilbert)*  makes  use  of  the  same 
type  of  computing  oscillator,  it  will  be  about  the  same  as  the  other  two  in 
this  respect  sc  long  as  the  orientation  is  well  away  from  the  singular  point. 
The  dependence  on  error  increases  as  the  singular  point  is  approached, 
however,  until  multiplier  errors  dominate  the  solution. 

The  fact  that  the-  multipliers  oscillate  only  half  as  fast  in  the  quaternion 
method  is  of  the  first  importance.  There  is  no  limit  to  the  speed  at  which 
the  servos  may  be  required  to  run  in  the  Euler  angle  method,  though  Howe 
and  Gilbert*  show  surprisingly  good  results  when  servos  are  slewing. 

It  should  be  observed  that  the  bandpass  of  even  the  rather  old  servos 
used  in  this  simulation  is  so  good  that  it  contributes  nothing  to  the  errors, 
so  long  as  the  sei  vos  arc  operating  linearly.  The  real  problem  in  cervos 
is  rate  and  acceleration  limiting.  It  is  on  this  point,  then,  that  the  lower 
servo,  frequency  of  the  quaternion  method  shows  its  real  advantage. 


*  Howe  and  Gilbert.,  op.  cit. 
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As  to  the  amount  of  computing  equipment  required,  the  table  below 
will  serve  as  a  useful  means  of  comparison. 


-  — -  -  .  . 

’ 

. 

Method 

Basic  Equipment 

First 

Each 

Amplifiers 

Integrators 

Multiplications 

Vector 

Additional 

Vector 

Quaternion 

36 

4 

26 

9  indep 

9  dep 

Dir.  Cosine 

31 

6 

36 

9  dep 

9  dep 

Euler  Angle 

_ ±  . 

6 

24 

12  dep 

12  dep 

The  figures  for  Euler  angles  were  taken  from  Howe  and  Gilbert  with  addition 
of  another  loop  to  compute  sin  4*  and  cos  »}»,  which  would  be  necessary  if 
complete  vector  tra;  sformations  were  to  be  made.  It  may  be  seen  that  in 
amount  of  equipment,  the  advantage  lies  with  Euler  angles,  with  the  quater¬ 
nion  method  second,  well  ahead  of  direction  cosines.  The  last  two  colums 
give  the  number  and  type  of  additional  multiplications  required  to  transfer 
the  first  vector,  and  each  additional  vector.  In  the  all-important  area  of 
multiplications,  it  is  seen  that  the  quaternion  method  is  nearly  as  good  as 
direction  cosines  and  will  be  better  if  a  large  number  of  vectors  are  to  be 
transformed. 

In  ease  of  programming  and  reliability  there  is  not  much  to  choose  between 
quaternions  and  direction  cosines,  except  that  the  latter  takes  somewhat  more 
equipment.  In  the  Euler  angle  system,  some  thought  must  be  given  to  keeping 
the  inevitable  division  circuit  stable  and  possibly  protecting  the  associated 
amplifiers,  but  this  is  not  serious.  The  high  card  of  Euler  angles  is  that  they 
are  so  easy  to  interpret.  A.fter  all,  a  prime  function  of  a  simulator  is  to  tell 
the  operator  or  engineer  what  the  simulated  system  is  doing,  and  for  ready 
interpretation,  there  is  nothing  like  the  Euler  angle  system.  In  the  first  place, 
most  people  in  the  aircraft.fi eld  know  what  "pitch",  "roll",  and  "bank"  mean 
whether  they  know  differential  equations  or  not.  Cockpit  presentations  of 
attitudes  are  given  in  terms  of  gimbal  angles,  which  are  nothing  but  Euler 
angles.  There  is  a  tendency  to  overrate  the  advantage,  however.  It  is  a 
matter  of  experience,  and  on<-  '-a  learn  to  interpret  a  transformation  matrix 
with  only  a  —  •>  *  of  effort.  Similarly  one  can  learn  to  interpret  the 

-  -  i  7 
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Finally,  then,  a  simple  statement  of  the  comparison  is  this:  Euler  ; 

angles  are  convenient  for  interpret  uion,  but  for  accuracy,  they  cannot 
compete  with  either  direction  cosines  or  quaternions.  On  all  significant 
bases  of  comparison,  the  quaternion  method  appears  superior  to  direction 
cosines.  Of  course,  it  should  be  kept  in  mind  that  we  arc  considering  only 
coordinate  transformations  capable  of  unrestricted  rotation  about  any  axis. 

It  might  not  be  out  of  place  at  this  point  to  interject  a  few  remarks  about 
the  general  use  of  three-dimensional  coordinate  conversions  of  the  type 
considered  here.  It  goes  almost  without  S3ying,  that  such  a  coordinate 
conversion  would  not  be  used  except  in  a  very  large  and  involved  simulation. 
Another  way  of  putting  it  is  that  in  order  to  make  this  type  of  coordinate 
conversion  useful,  other  parts  of  the  problem  must  receive  a  similarly 
general  and  unrestricted  treatment.  This  is  very  rarely  done.  The  reason 
is  that  the  amount  of  labor  involved  is  such  as  to  be  justifiable  by  only  the 
most  overwhelming  technical  reasons.  Computer  capacity  is  not  the  limiting 
factor.  There  are  many  computer  installations  in  the  country  whose  capacity 
is  equal  to  the  largest  simulations  yet  attempted.  The  problem  appears  to  lie 
in  the  tremendous  amount  of  painstaking  detail  involved  to  set  up  the  problem, 
check  it  out  and  keep  it  working.  The  only  reasonable  answer  to  this  is  a 
policy  of  programming  problems  in  larger  pieces.  This  is  comparable  to  use 
of  subroutines  on  a  digital  machine.  Without  going  into  the  matter  at  length, 
it  seems  clear  that  coordinate  conversion  lends  itself  to  this  technique  probably 
better  than  any  other  part  of  the  problem.  Consider  for  example,  the  quater¬ 
nion  method  of  Figure  Z,  The  inputs  to  the  coordinate  conversion  are  the  three 
voltages  P,  Q,  and  R,  and  the  outputs  are  the  nine  direction  cosines.  There 
are  only  twelve  gains  which  would  have  to  be  changed  from  one  problem  to  the 
next,  and  these  are  the  P,  Q,  and  R  product  inputs  to  the  four  integrators. 

These  serve  to  establish  the  maximum  allowable  rate  of  rotation  and  the  scale 
factor  on  the  inputs  P,  Q,  and  R.  Nothing  else  in  the  entire  circuit  would  have 
to  be  changed.  The  computer  operator  would  not  even  have  to  know  how  the 
coordinate  conversion  worked.  It  would  be  converted  into  that  "Mack  box"  of 
which  we  are  all  so  fond.  By  using  this  sort  of  technique  throughout  the  problem, 
the  amount  of  labor  involved  in  large  simulations  could  be  cut  by  the  factor  ten. 
This  would  simplify  things  other  than  the  setup  procedure.  The  engineer  would 
no  longer  have  to  spend  so  much  time  deciding  what  is  negligible.  He  could 
start  everything  in  and  find  what  is  negligible  by  throwing  it  out  and  seeing  if 
it  changes  the  result,  which,  of  course  is  in  the  best  mathematical  tradition. 
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APPENDIX  A  i 

ORTHOGONAL  TRANSFORMATIONS 


Fundamental  to  the  study  o£  rigid  body  motions  is  the  determination  oi 
how  many  degrees  of  freedom  it  has.  Putting  it  another  way,  the  problem 
is  to  determine  how  many  numbers  one  must  specify  in  order  to  describe 
the  orientation  of  the  body.  In  order  to  do  this,  it  will  also  be  necessary 
to  give  a  more  exact  definition  to  the  term  “rigid". 

Assume  that  a  body  is  composed  of  a  large  number  of  elementary 
particles.  If  the  distance  between  the  i  JLh.  particle  and  the  i  th  particle 
is  constant  of  Ml  particles  i  and  j,  then  the  bedy  is  said  to  be  rigid. 

If  all  the  N  partic.  ?s  were  independent  of  each  other,  it  would  require  3N 
coordinates  to  speciq  them  all.  (Three  cartesian  coordinates  are  required 
to  specify  ihe  position  of  a  point.)  The  particles  are  not  all  independent, 
however.  In  fact  the  Position  c-1  any  particle  in  the  body  may  be  specified  by 
the  distances  to  any  three  non-collinear  joints  in  the  body. 
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The  points  I,  2  and  3  in  Figure  8  have  been  chosen  at  random,  the  only 
condition  being  that  they  do  not  lie  along  the  same  line.  By  the  rigid 
body  condition  that  r  j  ,  and  ^  are  constant,  the  position  of  the  i  th 
particle  is  fixed  once  the  positions  of  the  particles  i,  it  follows  that  the 
position  of  every  particle  in  the  body  is  specified  once  the  three  point?  are 
specified.  In  other  words,  the  position  of  the  body  is  specified  by  the  positions 
of  these  three  points.  Specifying  three  points  would  require  nine  coordinates 
if  all  the  points  were  independent. 

There  are  three  conditions  to  be  fulfilled  by  these  coordinates,  however, 
namely  the  prescribed  values  of  r^,  rj3  and  r  iy  Thus  six  coordinates  are 
required  to  specify  the  position  of  the  rigid  body.  Another  way  ol  saying  th>s 
is  to  say  that  the  rigid  body  has  six  degrees  of  freedom.  These  are  frequently 
divided  into  two  groups  called  translational  and  rotational  degrees  of  freedom. 
The  three  coordinates  used  to  specify  the  orientation  of  some  point  in  the  body 
(say  the  point  1  in  Figure  8)  m  the  xyz  coordinate  system,  may  be  called  the 
translational  coordinates,  while  the  three  ccordi  nates  required  to  specify  the 
relative  orientation  of  the  other  two  points  could  be  called  the  rotational 
coordinates.  The  translational  coordinates,  then,  are  associated  with  the 
motion  of  the  body  as  a  whole,  while  the  rotational  coordinates  are  associated 
with  the  orientation  of  the-  bouv. 

2,  Orthogonal  Transformations 

Consider  a  vector  r  which  has  components  x,  y  and  z  in  the  XYZ 
coordinate  system.  If  the  unit  vectors  along  the  X,  Y  and  Z  axes  are 
called  i,  j  and  k,  then  it  is  possible  tu  write  r  as  follows: 


ix  +•  jy  +  kz 


/ 1  >  v 


Plow  assume  some  coordinate  system  X’Y’Z’  which  lias  the  same  origin  as 
the  XYZ  system  Out  an  arbitrary  rotation  with  respect  to  it.  The  components 
of  r  in  this  system  are  x'  y'  and  z'  and  the  unit  vectors  along  the  three 
axes  are  i',  j*  and  k'.  The  vector  r  may  also  be  written 


r  =  i'  x'  +  j  *  y '  +  k’z'. 


(123) 


The  problem  is  to  determine  the  components  x’,  y'  and  z'  in  terms  of  x,  y  and  ? 


vV A  DC  TR  58  -3  7 


57 


and  the  relative  orientation  of  the  two  coordinate  systems.  This  process 
is  called  an  orthogonal  transformation. 

It  is  possible  to  write  the  unit  vector  i'  in  terms  of  its  components  in 
the  XYZ  system, 

i'  =  {i '  -  T5T  +  (i'*j)j  +  (i'.k)k.  (124) 

Since  all  these  vectors  have  unit  magnitude,  the  dot  product  of  two  is  simply 
the  cosine  of  the  angle  between  them. 

i ’  *  i  =  cos  [_  i1*  i  =  ai  l* 


l'.j  =  cos 


IS'' T  *  ai2« 


(125) 


T'-k  =  cos  /?-k  =  a13. 

The  same  process  may  he  applied  in  obtaining  j'  and  k'. 

y  =  +  (y  +7)f  +  (T’*k)k» 

k'  =  (k'-i5T  +  (k'*T5T  +  (k'.kik', 

so  the  entire  set  of  relationships  may  be  written: 

V  =  a  nT  +  a12T  +  a13k  , 


j'  =  a21*  +  a22j  +  a23k  ’ 


(126) 


k-  =  a31i  +  a32j  +  a33k  . 

It  is  possible  to  apply  an  exactly  similar  process  in  expressing  the  unit  vectors 
i,  j  and  k  in  terms  of  their  components  in  the  X'Y'Z'  system. 


i  =  ani'  +  a21j'  +  a31k'  , 


j  -  a12r  f  a22j'  +  a3£k  , 
k  =  +  a Z3y  +  a33k'  . 

Figure  9  shews  the  two  coordinate  systems  and  the  unit  vectors. 


(127) 
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Figure  9. 

It  is  now  possible  to  determine  the  components  of  the  vector  r  in  the  X'Y'Z' 
coordinate  system. 

x-  =  r  •  f1  =  a,,x  +  a^y  +  a y.z  , 

=  r  -T  "  +  az,y  +  *,>3*  ,  (128) 

Z'  =  7.  k'  =  a31x  +  a32y  +  . 

The  nine  quantities  ^  -  -  a33  are  called  the  direction  cosines.  They  provide 
the  means  of  transforming  a  vector  from  one  coordinate  system  to  another  and 
therefore  they  specify  the  orientation  of  the  X'Y'Z1  system  with  respect  to  the 
XYZ  system.  It  was  developed  earlier  that  only  three  parameters  were  necessary 
vo  specify  the  orientation  of  a  rigid  body.  Therefore  there  must  be  six  equations 
relating  the  direction  cosines  to  each  other.  It  will  be  noted  that  regardless  ' 
of  what  rotation  is  applied  to  the  coordinate  system,  the  length  of  any  vector 
must  remain  unchanged.  This  means  that 

{x')Z  +  (y1)2  +  (z1)2  =  x2  +  y2  +  z2.  (129) 


vVADC  TR  58-17 


59 


60 


Matrices  are  a  type  of  mathematical  entity  which  may  be  conveniently 
applied  to  the  problem  of  rigid  body  rotations.  The  rules  for  manipulating 
these  quantities  will  now  be  reviewed. 


3.  Properties  of  Matrices 

The  multiplication  of  a  matrix  by  a  vector  is  the  first  operation  of 
interest.  Symbolically,  this  is  represented  by 


=  (A)i 


(134) 


For  convenience,  the  x,  y  and  z  components  of  r  are  denoted  by  x^  and 
Xy  Note  that  a  vector  "r  may  be  viewed  as  a  matrix  of  only  one  column. 
The  equation  might  be  written  i 

I 


The  rule  for  performing  this  operation  is 


3  ‘ 


.3  =  1 


(135) 


(136) 


If  these  operations  are  carried  out,  a  set  of  three  equations  is  obtained  which 
is  identical  with  the  set  of  Equations  (128).  This  means  that  multiplication  of 
a  vector  by  matrix  using  the  multiplication  rule  above  represents  a  transfor¬ 
mation  of  that  vector  from  one  coordinate  system  to  another.  For  this  reason 
the  matrix  (A)  may  be  called  the  transformation  matrix. 

The  case  of  two  successive  rotations  is  an  important  one.  Let  the  first 
rotation-be-represented  by  a.  matrix  («B).-  T-hen  the  components  of  a-vector 
after  this  rotation  will  be  given  by  .  , 

V  *  I  Vj  •  (U7> 

j 
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If  the  second  rotation  is  represented  by  the  matrix  (A),  then  thi  components 


of  the  vector  after  this  second  rotation  would  be 


*i'  ’  -  l  *ik*k  * 

k 

Substituting  (137)  into  (138)  gives 


I  aiV  X  Vi  • 


(138) 


I 

i  'k 


Note  that  this  can  be  out  in  the  form  of  Equation  136. 


c."  =  S  C..X.  , 

1  L  1j  i 


(139) 


where 


c  .  =  )  a.,  b,  .  . 

iJ  L  kj 

k 


(140) 


Thus  the  two  rotations  may  be  r  eplaced  by  a  single  rotation  (C),  the  element 
of  which  may  be  computed  from  (140).  Symbolically, 


(C)  =  (A)  (B)  . 

It  can  be  seen  by  the  rule  of  Equation  (140)  that 


(141) 


(A)  (B)  *  (B)  (A)  , 

so  the  process  of  matrix  multiplication  is  not  ccmmutative.  The  process  of 
matrix  multiplication  is  associative. 


(A)  [(B)  (C)  ]  =  [(A)  (B/  ]  (C). 

The  matrix  (A)  was  used  to  transform  the  vector  r  into  the  vector  r'. 

It  is  of  interest  now  to  investigate  the  properties  of  the  matrix  (A) 
which  transforms,  r'  into  r.  The  elemc-iits  of  this  inverse  matrsx  are  desig 
1  ■  .  a  The  inverse  matrix  is  defined  by  the  following  equation. 


(A)*1  (A)r  =  r. 


(142) 


A' ADC  Til 


Doing  the  first  operation,  the  result  is 


=  )  a..x.  . 

1  L  D  J 


(143) 


Now  applying  the  inverse  transformation  to  this  gives 


r  -  \  v  v  ■ 


- 1  v  1  %  h  • 

i  j 

V  -  1  (laH  xj ' 


J 


Now  according  to  the  requirement  that  this  must  give  back  the  original  vector, 
x^1  1  =  Xp,.  This  will  be  true  only  if 


>.  v  aij  =  6jk- 


(144) 


/ 1  \  - 1 


This  shows  that  the  product  of  the  two  matrices  (A)  and  (A)'  will  be 


/‘  0 


(A)-1  (A)  =  0  1  0  =  (I) 


0  0  1 


(145) 


This  matrix  (l)  is  called  the  identity  matrix.  It  may  oe  easily  shown  from  the 
rules  of  matrix  multiplication  that  for  any  matrix  (Q)„ 


(I)  (Q)  =  (Q)  (0  =  (O) 


1 14b; 


Now  siziee  (A) ~  corresponds'!©  some  physical  rotation;  there  must  exist  some 
matrix  (R)  which  i"  the  inverse  of  (A)~^.  In  other  words,  there  must  be  an 
(R)  such  that 


(R)  (A)-1  =  (1). 

Mow  if  (R)  is  appl  ed  to  both  sides  of  Equation  (145),  the  result  is 


(147) 


(R)  (A)-1  (A)  =  (R)  (I)  ^  (R). 


(148) 


WADC  TR  58-17 


63 


Since  matrix  multiplication  is  associative,  Equat.on  (147)  may  be  substitited 
into  Equation  (148)  to  give 


(1)  (A)  =  (R). 

(A)  =  (R) . 

This  means  that 

(A)'1  (A)  =  (A)  (A)'1  =  (I)  ,  049) 

so  that  (A)  and  (A)"*  commute. 

Now  consider  the  double  sum. 


1  ‘ijV  * 

k.i 


(150) 


This  sum  may  be  written  two  ways,  depending  on  the  order  of  summation. 


(151) 


Applying  Equation  (144)  to  the  quantity  in  parentheses  on  the  left  hand  side, 
and  applying  the  orthogonality  condition  of  Eqxiation  (131)  to  the  quantity  in 
parentheses  on  the  right  hand  side,  the  result  is 


I  6jk  V  =  I  6ii  aij  ’ 


(152) 


This  is  the  important  result.  To  form  the  inverse  of  an  orthogonal  matrix, 
the  rows  and  columns  are  simply  interchanged.  Note  that  this  conclusion 
holds  true  only  for  orthogonal  matrices.  This  is  because  the  orthogonality 
conditions  were  used  to  prove  Equation  (152).  In  general,  the  matrix  formed 
by  interchanging  rows  and  columns  is  called  the  transposed  matrix  and  is  desig 
nated  by  (A).  The  complex  conjugate  of  this  transposed  matrix  is  called  the 

-f  ^  ♦ 

adjoint  matrix  and  is  indicated  by  (A)  =  (A)  .  A  matrix  is  said  to  be  unitary 

if  it  satisfies  the  condition, 

(A)+(A)=(l).  (153) 

Of  course  tht.c-e  la  :--r  definitions  are  relatively  meaningless  in  the  case  of 
real  m;  cs.  However,  use  is  sometimes  made  of  matrices,  the  elements 
r-f  .  are  complex  numbers. 
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It  is  of  interest  to  investigate  the  characteristics  of  the  determinant 
formed  by  the  element  of  a  matrix.  T  to  determinant  of  the  matrix  (A) 
will  be  written  as  [A].  It  will  be  noted  that  the  law  of  matrix  multiplication 
is  the  same  as  the  law  for  multiplication  of  determinants.  Therefore, 

[ABj  =  [A]  [B].  (154) 

Evidently  the  determinat'd  of  the  identity  matrix  has  the  value  unity, 
therefore,  from  Equation  (145)  it  may  be  seen  that 

[A'1]  [A]  =  1,  (155) 

provided  that  (A)  is  orthogonal.  Since  interchanging  rows  and  columns  does 
not  alter  the  value  of  a  determinant,  [A-*]  =  [A]  and,  from  Equation  (155), 

[A]2  =  1.  (156) 

This  means  that  the  determinant  of  the  transformation  matrix  can  have  only 
the  values  plus  or  minus  one.  If  the  rotation  is  a  real  one,  it  may  be  shown 
that  +1  is  the  only  allowable  value.  There  is  a  certain  type  of  matrix  opera¬ 
tion  which  is  called  a  similarity  transformation.  It  is  defined  by 

(A)'  =  (B)  (A)  (B)'1.  (157) 

It  can  easily  be  shown  that  the  determinant  of  (A)  is  the  same  as  the  determinant 
of  (A)',  that  is  to  say,  the  value  of  the  determinant  of  a  matrix  is  invariant  under 
a  similarity  transformation  of  that  matrix.  This  may  be  shown  by  simply  apply¬ 
ing  both  sides  of  (157)  to  the  matrix  (B). 

(A)'  (B)  =  (B)  (A)  (B)"1  (B)  =  (B)  (A).  (158) 

From  this  it  is  seen  that 

■  ■  [A'J  fBl  -  I®]  (AJ  - . 059)  ' 

Since  [B]is  a  number  and  not  zero,  it  is  possible  to  divide  both  sides  by  it 
and  obtain  the  result 

[A]  =  [A'J  .  (160) 

which  demonstrates  the  proposition. 
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There  is  another  set  of  relationships  among  the  direction  cosines  which 
will  prove  to  be  of  interest.  Consider  the  set  of  Equations  (126).  If  the 
i',  j'  and  k'  vectors  are  mutually  perpendicular,  then  the  following  relation¬ 
ships  apply: 

i1  x  j'  =  k'  ,  j'  x  k'  =  i"  ,  k'  x  ?  =  j'  .  (161) 


If  these  vector  equations  are  expanded  in  the  unprimed  system,  and  their 
components  equated,  the  result  is 


all 

— 

a22a33"a23a32  * 

a21 

= 

al 3a32”a12a33  * 

a31 

= 

a12a23‘a13a22  * 

a12 

= 

a23a31 “al 1 a23  » 

a22 

= 

al la33"a13a31  * 

a32 

= 

al 3a21 _al 1 a23  » 

al  3 

= 

a21a32"a22a31  * 

a23 

= 

a12a31_alla32  * 

a33 

aila22"a12a21  * 

(162) 


These  nine  equations  are  really  consequences  of  the  orthogonality  conditions. 
They  present  a  means  for  solving  for  any  direction  cosine  in  terms  of  the  others. 


4.  Infinitesmal  Rotations 

t 

It  would  be  a  great  advantage  if  a  vector  could  be  associated  with  a  finite 
rotation,  but  it  turns  out  that  this  is  not  possible.  For  one  thing,  finite  rotations 
are  not  commutative;  nor  even  anti-commutative.  That  is  to  say  the  order  of 
the  operations  must  be  preserved.  While  this  is  true  of  a  finite  rotation,  it  will 
be  shown  that  a  vector  may  be  associated  with  an  infinitesmal  rotation  and  that 
therefore,  the  known  characteristics  of  vectors  may  be  used  in  the  treatment  of 
such  rotations.  Consider  the  matrix  that  describes  a  rotation  thru  the  angle 
'>"•  •*-  t'  *  ■  >-hich  makes  the  angles  a^,  Pj  and  Yj  with  the  X,  Y  and  Z 
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axes  respective1^',  This  matrix  may  be  potten  by  substituting 


^1  ^1 


J-  ,  cos  ~Y~  =  0 


into  the  matrix  (12)  and  dropping  higher  order  terms.  The  result  is 


1  Ap.j  cos  Yj  -  Ap.j  cos  [ij 


(A)j  =  I  -  AjJ^  cos  1 


Ap^  cos 


Ap cos  pj  -  Ap-j  cos  a.j  I 


(163) 


This  matrix  differs  only  slightly  from  the  identity  matrix.  This  may  be 
seen  more  clearly  by  writing  it  in  the  following  form: 


1  0  0 


(A)^  =  0  1  0  +  Apj  cos  Yj 


0  0  1/ 


Af-!  C°R  Yj 


Afij  cos  Pj  -  Au|  cos 


Afij  cos  Pj' 


Ap^  COS  Qj 


(164) 


This  latter  matrix  is  anti -symmetric  or  skew -symmetric.  Notice  that  this 
matrix  has  only  three  independent  elements,  Apj  cos  ;  Ap^  cos  Pj  ; 

Apj  cos  Yj  and  that  they  are  simply  the  three  components  of  a  vector  of 
magnitude  Ap  which  is  oriented  along  the  axis  of  rotation.  It  will  be  shown 
that  this  is  the  vector  which  may  be  associated  with  imintosrnal  rotation.  Let 
t’nesi  three  components  be  called  so  that  (101)  may  be  written 


(A)j  =  (I)  +  -  «3 


(165) 


No'A'  if  the  infinitesmal  rotation  (A)^  is  followed  by  another  infinitesimal  rotation 


(a)’j  ,  of  ilie  form 


(aJ'j  =  (i)  *  -  n’3  o 


-  n\ 


(166) 
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tiun  the  combined  relation  (A)^  (A)j  5s  seen  to  be  the  following,  if  higher 
unlet  infinitesimals  are  dropped; 


(A)''  =  (A)’  (A)  =  (I)  +  -  S2-  0 


0  PJ' 


Q"?  -  V'l 


(167) 


where 


Q1^  =  +  Q3  ;  =  Piy  +  Oa  ;  =  Q'{  +  flj  . 

Since  the  second  order  infinitesmals  were  dropped,  the  order  or  sequence  of 
the  infinitesmal  rotations  is  unimportant.  This  is  one  condition  which  is  neces 
sary  if  these  rotations  are  to  be  represented  by  vectors.  From  the  makeup  of 
,  fi'L  ,  it  is  seen  that  the  vector  representing  the  combined  rotation  is 
simply  the  sum  of  the  two  vectors  tor  the  single  rotations, 

A  more  conclusive  demonstration  of  the  fact  that  the  quantities  S7j,  ft 

are  the  components  of  a  vector  associated  with  the  infinitesmal  transformation 
is  the  demonstration  that  the  matrix  components  transform  like  components  of 
a  vector  under  a  coordinate  transformation.  Consider  ,i  matrix  (A)  which  oner 
ates  on  a  vector  R  to  produce  a  vector  R', 


R'  =  (A)R. 


(168) 


Now  if  an  additional  matrix  (B)  is  applied  to  this  equation, 

(B)H'  =  (B) (A)R, 

(B)R'  =  (B)  (A)  (FI)-1  (>3)R. 


(169) 


This  equation  is  simply  Equation  (168)  when  seen  iu  a  different  coordinate 
system,  and  v  o)  (^w  (13)  is  tl'iC  matrix  (y\)  when  viewed  from  the  different 
coordinate  system.  This  is  the  similarity  transformation,  which  has  been 
introduced  before.  If  a  similarity  transformation  is  applied  to  the  'matrix  of 
Equation  (16  5),  the  result  is 


(15)  (A) j  (131  “  =  (A)’  -  (I)  t  ..  0 


0  fi*. 


(170) 
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Expanding  and  equating  components 


*  bll 

+ 

b12 

J>2 

+ 

b13 

«3. 

* 

^  b2i 

"l 

+ 

b22 

°2 

i 

b23 

S73, 

^3 

=  b31 

"l 

+ 

b32 

*2 

h 

b3  3 

(171) 


Thus,  the  infinitesmal  transformation,  when  viewed  from  the  other  coordinate 
system  defined  by  (D)  is  still  nearly  the  identity  transformation,  and  the  vector 
which  represents  the  vector  associated  with  the  infinesmal  transformation  in 
this  new  system  is  simply  the  transform  of  the  vector  representing  the  infini¬ 
tesimal  transformation  in  the  other  coordinate  system.  'I  his  shows  the  vector 
character  of  the  set  of  elements  Cl j, 

By  using  this  infinitesmal  transformation,  the  rate  of  change  of  the  trans¬ 
formation  matrix  (A)  may  be  found  in  much  the  same  way  that  the  derivative  of 
the  matrix  (H)  was  established  in  Section  J.IL  If  (A)  is  the  matrix  at  the  begin¬ 
ning  of  time  interval,  and  (A)  is  the  matrix  at  the  end  of  time  At,  then  the 
derivative  of  (A)  is  given  by 


-jr  (A)  -•  lim 

dt  A-0 


(A)1  -  (A) 
At 


]• 


(172) 


(A)1  may  be  viewed  as  the  rotation  (A)  followed  by  the  infinitesmal  transform¬ 
ation  going  from  (A)  to  (A)'.  In  other  words 


where 


(A)'  =  [(I)  +  (t)]  (A), 


0  cos  y  -  cos  p1 
(«)  =  Ap  (-  cos  V  0  cos  a 

i  cos  p  -  cos  a  0  i 

Thus  the  derivative  becomes 


dl  (A)  =  lim-j^(A)  liraj- 
ai  A-0  **  A-0  ^ 


COS  p 


cos  y 
0 

-  cos  o 


(173) 


(174) 


Again,  in  the  limit is  simply  the  rate  of  rotation,  and  cos  o  =  P, 


At 

-gj~-  cos  P  =  Q,  cos 


»  y 


R,  so  the  equation  becomes 
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a32 

a33  /  \Q 

-P 

\*31 

a32 

ai3 

Expanding,  and  equating  components  gives 


all 

5 

a21R 

- 

a31 

a12 
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a22R 

- 

a32 

a13 
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a23R 

- 

a33 

a21 

= 

a31P 

- 

all 

a22 

= 

a32P 

- 

al  2 

= 

a33P 

- 

a13 

a31 

= 

allQ 

- 

a21 

a32 

s 

a12° 

- 

a22 

• 

a33 

= 

a13A 

- 

a23 

(176) 


These  are  the  rates  of  change  of  the  direction  cosines  in  terms  of  the  angular 
velocity.  Now  if  Equation  (175)  be  multiplied  on  the  right  by  the  transpose  of 
(A),  the  result  is 


/ail 

a12 

'l>\ 

/a, . 

/  11 

a21 

a31^ 

(° 

R 

a21 

a22 

[  a12 

a22 

a32 

= 

-R 

0 

'31 

a32 

W 

\13 

a23 

a33^ 

-P 

(177) 


Expanding  and  equating  components  gives  the  following  relationships: 


P 

— 

a31 

a21 

+ 

a32 

a22 

+ 

a33 

• 

a23  ’ 

_p 

- 

a21 

a31 

+ 

a22 

a32 

a23 

• 

a33  • 
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= 
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a31 

+ 

a12 

a32 

+ 
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a33  » 

-Q 

a31 
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al  3  ’ 
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a21- 

kll 

+ 

a22 

a12 

+ 

a23 

• 

al  3  » 

-R 

= 

all 

a21 

+ 

a12 

a22 

+ 

a13 

• 

a23  * 

(178) 
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It  is  interesting  that  two  different  expressions  are  obtained  for  each  of  the 
velocity  components.  This  is  a  consequence  of  the  great  amount  of  redundancy 
in  the  direction  cosines.  The  equivalence  of  the  two  expressions  for  any  one  of 
the  components  may  be  shown  by  making  use  of  Equation  (162). 
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APPENDIX  B 


THE  EULER  ANGLES 

It  was  demonstrated  in  Appendix  A  that  three  parame'.ers  were  required 
to  fix  the  orientation  of  a  rigid  body  and  hence  of  a  coordinate  system.  The  nine 
direction  cosines  do  not  lend  themselves  to  a  reduction  to  three  simple  para¬ 
meters,  nor  do  they  give  a  very  lively  picture  of  the  orientation  of  the  body. 

Both  these  difficulties  are  overcome  by  use  of  Euler  angles,  the  only  three-para¬ 
meter  system  in  common  use.  In  this  method,  a  rotation  is  represented  by  three 
individual  rotations  taken  in  a  specified  sequence  about  certain  specific  axes.  In 
the  literature,  there  is  no  agreement  whatever  on  the  order  of  rotations,  the  axes 
about  which  the  rotations  are  made,  or  notation.  These  are  varied  to  suit  the 
needs  of  the  problem  and/or  the  author's  whim.  Texts  on  classical  mechanics 
give  sets  of  angles  defined  so  as  to  facilitate  solution  of  the  spinning  top  problem. 
The  system  presented  here  is  the  most  common,  though  by  no  means  the  only  one 
used  in  aircraft  work. 

Consider  two  coordinate  systems  initially  coincident.  One  set  of  coordinates, 
the  x,  y,  z ,  will  be  referred  to  as  the  fixed  system,  and  the  other  will  move  with 
respect  to  it.  The  first  rotation  is  through  the  angle  ^  about  the  z'  axis.  This 
is  shown  in  Figure  10. 
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THIRD  ROTATION  -  ROLL  ANGLE  $ 


Figure  \d 


--n 


>:•  n.  ir<  gr 


The  second  rotation  is  through  the  angle  0  and  is  done  ibout  the  Y’  axis  and 
•  lie  resulting  axis  system  is  called  X",  Y",  Z".  This  rotation  is  shown  in 
Figure  11.  0  is  commonly  called  the  pitch  angle.  'lie  final  rotation  is  done 

about  the  X"  axis  through  the  angle  .  This  is  called  the  roll  angle  and  all 
three  rotations  are  shown  in  Figure  12.  Note  that  all  three  of  these  rotations 
are  in  the  positive  sense.  That  is  to  say  if  the  thumb  of  the  right  hand  is  placed 
along  the  axis  of  rotation,  the  direc.tiun  of  rotation  is  that  direction  in  which  the 
curled  fingers  point. 

It  is  now  necessary  to  determine  the  transformation  matrix  in  terms  of  these 
Euler  angles.  It  was  shown  earlier  that  successive  rotations  could  be  represented 
by  a  matrix  which  is  a  product  of  the  matrices  of  the  individual  rotations.  It  ia 
necessary  then,  only  to  compute  the  matrix  corresponding  to  each  of  the  Euler 
angle  rotations  and  to  multiply  them  together  in  the  appropriate  order.  Note  that 
each  of  the  rotations  is  simply  a  two-dimensional  transformation  because  in  each 
case  the  rotation  is  about  one  of  the  moving  axes  and  hence  components  along  that 
axis  are  unchanged. 

Consider  first,  the  rotation  through  the  angle  which  is  shown  in  Figure  10„ 

If  this  is  viewed  from  ubo-  e,  the  transformation  of  some  arbitrary  vector  R  would 
appear  as  shown  in  Figure  13 


/  *  ' 
I  1 


r  r 

L  7 
' . A 


V  \L>C  TR  F  8  - 1  7 


I  igure  13 


<1?9) 


It  can  be  seen  from  the  geometry  of  Figure  13  tli.it  the  new  x*  end  y1 
components  .ire  related  to  the  old  by  the  e<  nations 

x'  =  x  cos  41  +  y  sin  4*  , 
y'  =  -x  sin  3  y  cos  4‘  . 

Since  the  rotation  was  about  the  z  axis,  any  z  component  of  R  would  remain 
unchanged.  In  other  words,  Z  -  Z',  This  fact,  together  with  the  Equation  (179) 
shows  that  the  matrix  for  the  rotation  is 


/  cos  4*  sin  4  0\ 

1-sin  4*  cos  4'  0 

0  0  li 


(180) 


Now  the  rotation  of  Figure  11  may  be  viewed  from  the  front  along  the  Y’  axis, 
and  Figure  14  is  obtained. 


Figure  14 

From  the  gee  met  ry  of  the  above  figure,  it  may  be  seen  that 


-x1  cos  0  -  z'  sin  0, 
-x’  sin  0  3-  z’  cos  6  • 


(181) 
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In  tins  rotation,  the  Y  components  remain  unchanged  so  that  Y"  =  Y'. 
Tliereloro,  the  matrix  for  this  rotation  is 


c  o  s  0  0  -sin.. 

0  1  0 
sin  0  0  cos  0 


(18?.) 


The  final  rotation  may  be  viewed  from  the  front,  looking  along  the  X"  axis  of 
Figure  12. 


I 


From  the  geometry  of  this  figure  it  is  seen  that 


-Ml  - 


y"  COS  <f>  -t  z"  sin  A, 
-y"  sin  4*  +  z:' '  cos  if. 


(183) 


In  this  rotation,  the  X  components  remain  unchanged  so  that  X'" 
Therefore,  the  transformation  matrix  for  this  rotation  is 


X' 


0 

cos  i 
-sin 


0 

sin  y 

cos  / 


(184) 
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In  order  to  get  the  total  transformation  matrix  which  results  from  these  three 
rotations,  it  is  only  necessary  to  multiply  Me  three  individual  matrices  in  the 
correct  order. 


h 

0. 

°  ) 

1  cos  0 

0 

-sin  0 \ 

|cos  4j 

sin  4* 

O' 

0 

cos  <}> 

sin  <|> 

0 

L 

0 

-sin  4 

cos  4* 

0 

U 

-sin  <j> 

cos  4>/ 

\  sin  0 

0 

cos  0  | 

l  o 

0 

1, 

cos  0  cos  v(i 
-  sin  4*  cos  4> 

+  sin  <}>  sin  0  cos 


cos  0  sin  4< 
cos  4*  cos  4* 

+  sin  4>  sin  0  sin  4* 


-sin  0 


sin  <}>  cos  0 


.  +  sin  4>  sin  4> 

\,+  cos  4>  sin  0  cos  4* 


-sin  4>  cos  4< 

+  cos  <{>  sin  0  sin  ^ 


cos  0  cos  <)> 


(185) 


By  comparison  of  this  matrix  with  the  matrix  (A)  it  may  be  seen  that  all  of  the 
direction  cosines  and  hence  the  complete  transformation,  can  be  expressed  in 
terms  of  the  three  independent  parameters^  4**  9»  <K 

Since  the  position  of  a  coordinate  system  may  be  specified  in  terms  of  Euler 
angles,  the  rate  of  rotation  of  that  coordinate  system  must  be  related  to  the  rates 
of  change  of  the  Euler  angles.  We  now  investigate  this  relationship. 

It  is  shown  in  Appendix  A  that  a  vector  could  be  associated  with  a  rate  of 

rotation.  This  vector  is  along  the  instantaneous  axis  of  rotation  and  is  equal  in 

magnitude  to  the  rate  of  rotation.  Thus,  each  of  the  Euler  angle  rates  may  be 

associated  with  a  vector  along  the  axis  of  rotation.  Observe  that  the  vector 

associated  with  the  4*  rotation  of  Figure  10  is  directed  along  the  Z  axis  and 

points  downward  if  4*  is  positive.  Similarly,  the  rate  of  rotation  due  to  the  0 

rotation  of  Figure  11  is  a  vector  along  the  Y1  axis,  and  if  0  is  increasing,  the 

vector  is  in  the  positive  y'  direction.  Finally,  a  positive  roll  rotation  is  a  vector 

-directed  along  the  X'"  axis-of  Figure  12.  The  three  vectors  .representing  the  three 

individual  Euler  angles  rates  must  be  added  together  in  order  to  get  the  entire  rate 

of  rotation  of  the  system.  Recall  that  these  vectors  are  *.dded  according  to  the  usual 

vector  rule.  The  situation  is  shown  in  Figure  16  where  all  the  Euler  angle  rates  are 

assumed  positive.  Note  that  these  three  vectors  are  not  mutually  orthogonal.  The 
*  •  •  a 
4*  vector  is  normal  to  the  0  vector,  and  the  0  vector  is  normal  to  the  4>  vector, 

•  •  * 
but  the  <}>  vector  is  not  normal  to  the  4*  vector.  In  any  case,  the  three  may  be 
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transformed  into  the  X1"  YMI  Z1"  and  added  to  give  the  entire  velocity  vector, 
the  i)j  vector  has  the  components  0,  0,  4*  in  the  XY Z  system,  so  to  transform 

this  into  the  X'"  Yl,:  Z"1  system,  it  is  necessary  to  apply  the  full  transformation 
matrix  (185)  to  this  vector.  If  this  is  done,  the  result  is 

vj/'"  =  -i"'  4;  sin  0  +  j1"  4>  sin  4>  cos  0  +  k1”  4>  cos  0  cos  4>  *  (186) 

Now  the  vector  0  has  the  components  0,  0,  0  in  the  X"  Y"  Z"  coordinate  system. 
In  order  to  get  this  into  the  X,M  Y1"  ZIM  system,  it  is  only  necessary  to  transform 
through  the  last  of  the  Euler  angle  rotations  which  is  defined  by  the  matrix  (184).  If 
this  is  done,  the  result  is 

0"'  =  7"'  0  cos  4>  -  K'"  9  sin  4>.  (187) 

The  vector  4>,  of  course,  is  already  in  the  X"'  Y'”  Z"'  system,  being  defined  by 

4>  ^  7'"  4>  .  (188) 


In  order  to  get  the  entire  velocity  vector,  it  is  only  necessary  to  add  the  last  three 
equations.  If  this  is  done,  and  if  the  total  angular  velocity  vector  is  defined  as 

7  =  7"  P  +  7'"  Q  +  K"‘  R,  then, 

•  * 

P  =  4»  -  4»  sin  8, 

Q  =  0  cos  +  4*  sin  cos  0,  (189) 

•  • 

R  =  4^  cos  9  cos  4>  -  0  sin4>. 

•  •  • 

These  three  equations  may  be  solved  for  4*»  9.  41  giving 


.V  _  r,  cos  4)  ,  sin  <*> 

*  “  R  75T7  +  Q  7570  * 


9  =  Q  cos  <}>  -  R  sin  '9  ,  ~ 

4>  =  P  +  Q  tan  0  sin  <j>  +  R  tan  9  cos  4> 


(190) 


From  these  equations,  it  is  easier  t--'  see  the  difficulties  which  arises,  when  9 

o  •  • 

approaches  90  .  For  t'-is  v-’_e  of  9,  both  4*  and  4>  are  infinite.  It  is  interesting 
to  no,r  '  '  -ve  of  0  itself  has  no  such  anomalies. 
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APPENDIX  C 

SERVOMULTIPLIER  ERRORS 

In  the  past,  it  has  evidently  not  been  the  practice  to  treat  analog  multiplier 
errors  with  much  care  in  error  analyses.  The  usual  procedure  has  been  to 
assume  some  upper  bound  for  the  error  and  to  consider  the  error  constant  at  that 
value.  This,  of  course,  gives  a  pessimistic  view  of  the  results  to  be  expected, 
though  in  many  cases  this  is  not  undesirable,  and  it  certainly  gives  an  estimate  of 
the  order  of  magnitude  of  the  resultant  error.  Possibly  the  principal  reason  that 
errors  are  not  treated  with  more  care  is  that  little  is  known  about  their  detailed 
character.  No  two  multipliers  are  the  same,  a  given  multiplier  will  change  with 
time,  and  it  is  a  great  bother  to  measure  the  errors  anyhow.  It  would  appear, 
then,  that  a  statistical  approach  is  indicated.  If  the  statistical  distribution  of 
errors  can  be  establi.hed,  then  it  can  be  used  to  predict  the  distribution  of  errors 
in  the  problem  result.  This  will  not  give  a  measure  of  the  error  in  any  particular 
run,  but  it  is  not  practical  to  do  this  in  any  event.  It  would  appear  that  if  we  know 
the  distribution  of  errors  in  the  solution,  it  is  fair  to  say  that  we  know  all  that  is 
necessary  about  those  errors. 

It  will  become  clear  that  the  statistical  procedures  applied  in  this  work  are  of 
.the  very  simplest  sort.  Many  more  things  could  have  been  dene,  even  with  the  data 
which  were  taken,  but  it  was  not  felt  to  be  worthwhile.  The  prinicpal  reason  for  this 
is  that,  as  was  pointed  out  earlier,  multiplier  errors  did  not  turn  out  to  be  the 
principal  source  of  drift  in  the  coordinate  conversion  simulation.  Therefore,  it  was 
not  possible  to  check  the  predicted  error  distribution  against  the  observed  distribu¬ 
tion,  even  to  the  extent  permitted  by  the  small  number  of  samples  available. 

Generally  speaking,  the  procedure  was  to  take  an  average  of  all  errors  of  all 
multipliers  over  their  entire  range;  taking  the  average  in  a  maimer  generally 
consistent  with  the  way  the  multipliers  are  used  in  the  coordinate  conversion  simu¬ 
lation.  It  was  assumed  that  the  error  voltage  was  directly  proportional  to  the  voltage 
across  the  multiplying  potentiometer,  and  this  was  kept  at  ±100  volts  throughout 
the  measurements.  Once  the  error  distribution  is  established  for  this  case,  it  is 
possible  to  get  the  distribution  for  any  other  pot  voltage  by  simply  multiplying  the 
±100  volt  distribution  by  V  / .00,  where  V  is  the  pot  voltage. 


! 
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T.'.o  type  of  measurements  made  can  perhaps  best  be  i) ] us tr 2. ted  by 
cor.s.ee'.itior.  of  the  schematic  of  Figure  17,  This  shows  the  circuit  used 
to  test  one  serve,  with  all  its  pots. 


i  C.r  -  95  V 


Figure  17 

±100  volts  are  applied  across  all  pots.  The  wipers  arc  swept  linearly  from  one 
end  to  the  other  by  use  of  the  integrator.  It  starts  at  -95  volts  and  goes  to  ±95 
volts  at  the  rate  of  approximately  5  voits/ second.  The  entire  sweep  requires 
about  40  seconds.  At  this  slow  speed,  the  servo  tracking  error  is  essentially  xero. 
It  is  approximately  three  times  smaller  than  the  pot  errors.  During  the  sweep,  the 
voltage  on  each  of  the  multiplying  pots  is  being  compared  with  that  on  the  feedback 
pot,  ar.d  the  difference  is  recorded.  All  pots  are  loaded  with  0.1  megohms.  It 
might  appear  somewhat  more  logical  to  compare  the  multiplier  pot  voltage  with  the 
input  rather  than  with  the  feedback  pot  voltage.  Comparison  with  the  input  would 
include  the  tracking  error.  Even  aside  from  the  fact  that  the  tracking  error  is 
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small,  it  is  actually  desirable  to  exclude  it.  The  point  at  issue  here  is  the 
static  accuracy  of  the  multipliers.  The  tracking  error  is  associated  with  the 
dynamic  performance,  which  is  treated  by  different  methods.  In  the  steady 
state,  there  will  be  little  or  no  tracking  error,  and  the  only  error  remaining 
will  be  the  difference  between  the  feedback  pot  and  the  multiplying  pot.  This 
is  exactly  what  is  being  measured  by  the  circuit  indicated. 

A  sample  record  taken  in  this  way  is  shown  in  Figure  18.  This  is  the  one 
for  Servo  1,  and  shows  many  features  typical  of  the  others  as  well.  Observe 
that,  especially  on  Fbt  A  a  periodicity  may  be  seen.  There  are  ten  cycles  of 
the  oscillation,  so  it  is  evidently  associated  with  the  ten  turns  of  the  pots  used 
in  the  multiplier.  Not  all  records  show  this,  but  many  do.  It  may  also  be  seen 
that  there  is  a  rough  correlation  of  all  four  traces,  that  is,  they  all  tend  to  have 
the  same  sign  at  the  same  point,  and  the  peaks  on  all  traces  are  nearly  coincident. 
This  is  probably  due  to  the  fact  that,  for  this  servo,  the  feedback  pot  is  more 
nonlinear  than  the  multipli'r  pots.  On  many  of  the  records  this  correlation  was 
not  observed,  indicating  tha.  the  feedback  pot  was  the  more  accurate.  Also,  it 
is  seen  that  there  are  some  poor  contact  areas  on  Pot  A,  though  this  was  en¬ 
countered  in  only  a  few  cases.  It  should  be  mentioned  that  records  of  the  type 
of  Figure  18  are  taken  for  all  servos  periodically  by  the  Analog  Computation  Branch, 
Aeronautical  Research  Laboratory.  The  principal  purpose  of  this  is  to  check  align¬ 
ment  of  the  pots,  and  check  for  noisy  spots  such  as  appear  on  Pot  A  of  Figure  18. 

In  this  connection  it  should  be  pointed  out  again  that  no  special  adjustment  of  the 
multipliers  was  made  for  this  simulation.  The  error  records  used,  of  which 
Figure  18  is  an  example,  were  not  made  especially  for  this  study,  but  were  made 
about  two  months  earlier,  as  part  of  normal  maintenance  of  the  equipment.  Thus 
the  results  may  be  considered  typical  of  what  may  be  obtained  on  this  type  of 
equipment  in  normal  use.  Especially  in  view  of  this,  the  multiplier  accuracy  does 
seem  to  be  quite  good. 

On  each  of  the  error  recordings,  approximately  45  equally  spaced  readings 
were  made.  For  each  reading,  the  value  of  the  error  trace  was  sorted  into  class 
intervals,  rather  than  being  read  as  a  number.  This  was  done  to  simplify  the 
reading  process,  since  it  was  felt  that  all  desired  results  could  be  obtained  by 
dealing  with  frequencies  in  class  intervals.  All  intervals  were  of  equal  width 
(0.  02  volts).  Table  3  gives  the  results  obtained.  The  error  voltage  indicated  at 
the  top  of  each  column  is  t^e  cor1.,  of  the  class  interval. 
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The  frequencies  for  all  pots  of  all  servos  were  combined  into  a  single 
distribution  with  equal  weight.  This  process  could  be  justified  by  either  of  two 
assumptions:  the  error  voltage  is  completely  uncorrelated  with  shaft  position, 
or  all  shaft  positions  are  equally  probable.  The  first  of  these  has  to  do  with  the 
nature  of  the  servo,  and  the  second  with  the  nature  of  the  problem  being  solved. 

It  would  appear  that  both  assumptions  hold  in  the  present  case.  As  to  the  first, 
the  error  appears  to  be  uncorrelated  with  shaft  position  so  long  as  the  padding 
resistors  are  correctly  adjusted.  From  the  records  taken,  it  appears  that  they 
were  so  adjusted  for  the  servos  used.  Any  remaining  correlation  would  have  to 
be  repeatable  characteristic  of  the  pot  winding  machine,  and  no  such  character¬ 
istic  appears  on  cursory  examination  of  the  records.  Of  course,  the  10-cycle 
periodicity  appears  on  some,  but  not  on  others,  so  it  was  concluded  that  there 
was  no  important  correlation  between  shaft  position  and  error  voltage.  Matters 
of  this  kind  should,  of  course,  be  examined  by  application  of  statistical  techniques 
to  the  data,  and  in  fact  a  start  was  made  in  this  direction.  For  the  reason  cited 
earlier,  however,  it  wa.  decided  that  the  simple  procedure  described  would  be 
adequate. 

As  regards  the  second  assumption  above,  fhat  all  shaft  positions  are  equally 
probable,  this  appears  to  be  true  in  a  general  *ay  for  the  coordinate  conversion 
simulation.  For  other  applications,  of  course,  the  situation  might  be  different. 

In  either  thf»  quaternion  or  direction  cosine  method  it  is  true  that,  regardless  of 
the  orientation,  some  of  the  quantities  will  be  large  and  others  small.  A  complete 
rotation  of  the  coordinate  system  causes  certain  of  the  multipliers  to  sweep  through 
their  entire  range.  From  such  considerations,  the  hypothesis  that  all  shaft  positions 
are  equally  likely  appears  reasonably  sound. 

All  the  data  ot  Table  3  were  combined  into  a  single  cumulative  distribution 
which,  is  shown  as  the  solid  curve  of  Figure  19.  This  has  been  plotted  on  cumula¬ 
tive  probability  paper  so  that  a  normal  distribution  would  appear  as  a  straight  line. 

It  may  be  seen  that  this  curve  does  not  pass  through  the  50  per  cent  point  at  zero 
error  (its  mean  is  not  zero)  and  that  it  departs  from  linearity  for  positive  errors. 
This  suggests  a  bias  in  the  data.  An  examination  of  the  data  of  Table  3  shows  that  * 
Servo  l  has  something  wrong  with  it.  All  the  pots  have  a  large  positive  bias,  and, 
especially  Pot  A  has  an  unusually  large  dispersion.  If  the  data  for  this  servo  are 
deleted,  the  dashed  curve  of  Figure  1°  is  obtained.  It  may  be  seen  that  now  the 
curve  passes  through  the  50  n--.  cent  point  at  zero  error,  and  it  is  much  more 
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nearly  a  straight  line.  There  is  still  some  deviation  from  linearity,  and  this 
is  no  doubt  due  to  bias  in  some  of  the  pots,  but  it  appears  that  the  errors  may 
be  represented  reasonably  well  by  a  normal  distribution.  The  standard  devia¬ 
tion  may  be  determined  from  the  value  at  which  the  curve  crosses  the  84.  15  per 
cent  line.  This  turns  out  to  be  0.  053  volts.. 

In  the  error  analysis,  then,  it  may  be  assumed  that  the  multiplier  errors  are  (' 

| 

normally  distributed,  with  a  standard  deviation  of  0.  053  volts,  when  the  voltage 
across  the  multiplying  pot  is  ±100  V.  It  would  be  proportionately  less  for  smaller 
pot  voltages. 
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